NASA TN 0-5427 



FORTRAN PROGRAM FOR CALCULATING 
TRANSONIC VELOCITIES ON A 
BLADE TO-BLADE STREAM SURFACE 
OF A TURBOMACHINE 


by Theodore Katsanis 

Lewis Research Center 
Cleveland, Ohio 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • SEPTEMBER 1969 




1 . 


4, 


7 . 


9 . 


2 . 


Report No. 


2. Government Accession No. 


3. 


NASA TND-5427 


FORTRAN PROGRAM FOR CALCULATING TRANSONIC 
VELOCITIES ON A BLADE-TO-BLADE STREAM 
SURFACE OF A TURBOMACHINE 


Author( s) 


Theodore Katsanis 

P er fo rm i n g Organizotion Name ond Address 

Lewis Research Center 

National Aeronautics and Space Administration 


Cleveland, Ohio 44135 


Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 


Washington, D.C. 


20546 


14. 


Recipient's Catalog No. 

Report Date 

September 1969 

Performing Organization Code 
Performing Organization Report No. 

E-5046 

Work Unit No. 

720-03 

Contract or Grant No. 

Type of Report and Period Covered 

Technical Note 


Sponsoring Agency Code 


15. 


Supplementary Notes 


16 . ^ bs i ^ a e c t ' hod hag been deve i opec i to obtain a transonic flow solution on a blade-to-blade sur- 
face of a turbomachine. The blade may be fixed or rotating. The flow may be axial, 
radial, or mixed. The results include velocities on the blade surface and throughout the 
passage. The transonic solution is obtained by a velocity-gradient method, using infor- 
mation obtained from a finite-difference stream-function solution at a reduced weight 
flow. 


17. KeyWords (Suggested by 


A 


uth ot(s )) 


IB. 


Distribution Statement 

Unclassified - unlimited 


19. Security Classif. {of this report) 

Unclassified 


20. Security Classif. (of this page) 

Unclassified 


21« No. of Pages 

96 


22. Price" 


$3.00 


*For sale by the Clearinghouse for Federal Scientific and Technical Information 

Springfield, Virginia 22151 





CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION 1 

SYMBOLS 2 

MATHEMATICAL ANALYSIS 4 

NUMERICAL EXAMPLES 10 

Axial Stator 10 

Radial Inflow Turbine 14 

DESCRIPTION OF INPUT AND OUTPUT 18 

Input 19 

Instructions for Preparing Input 22 

Units of measurement 23 

Blade and Stream channel geometry 23 

Loss of relative stagnation pressure 23 

Inlet and outlet flow angles 23 

Defining the mesh 24 

Overrelaxation factor 25 

Format for input data 25 

Incompressible flow 25 

Straight infinite cascade 25 

Axial flow with constant stream- channel thickness 26 

Output 26 

ERROR CONDITIONS 38 

PROGRAM PROCEDURE 41 

Subroutine PRE CAL 43 

Calculation of A, W le and W^ e 43 

Calculation of W and maximum values of mass flow parameter pW 43 

Calculation of w, u>, and A for reduced weight flow 43 

Calculate j3. n and 0 out 43 

Subroutine TANG 43 

Subroutine TVELCY 44 

Subroutine VELGRA 44 

Subroutine BLOCKD 45 


iii 



Main Dictionary 

Program Listing For Subroutines Using Main Dictionary 47 

Subroutine CON TIN 

Subroutine ENTGRL g 2 

APPENDIXES: 

A - DERIVATION OF THE VELOCITY- GRADIENT EQUATION 88 

B - DEFINING THE REDUCED WEIGHT FLOW PROBLEM 92 

REFERENCES Q , 


FORTRAN PROGRAM FOR CALCULATING TRANSONIC VELOCITIES ON A 
BLADE-TO-BLADE STREAM SURFACE OF A TURBOMACHINE 
by Theodore Katsanis 
Lewis Research Center 

SUMMARY 


A method has been developed to obtain a transonic flow solution on a blade -to -blade 
surface between blades of a turbomachine. A FORTRAN IV computer program has been 
written based on this method. The flow must be essentially subsonic, but there may be 
locally supersonic flow. The solution is two-dimensional, isentropic, and shock free. 
The blades may be fixed or rotating. The flow may be axial, radial, or mixed, and there 
may be a change in stream channel thickness in the through-flow direction. A loss in 
relative stagnation pressure may be accounted for. 

The program input consists of blade and stream-channel geometry, stagnation flow 
conditions, inlet and outlet flow angles, and blade -to -blade stream -channel weight flow. 
The output includes blade surface velocities, velocity magnitude and direction at all in- 
terior mesh points in the blade -to -blade passage, and streamline coordinates throughout 
the passage. 

The transonic solution is obtained by a combination of a finite -difference, stream- 
function solution and a velocity- gradient solution. The finite-difference solution at a 
reduced weight flow provides information needed to obtain a velocity-gradient solution. 

This report includes the FORTRAN IV computer program with an explanation of the 
equations involved, the method of solution, and the calculational procedure. Numerical 
examples are included to illustrate the use of the program, and to show the results which 
are obtained. 


INTRODUCTION 

Two useful techniques for calculating blade surface velocities are the velocity- 
gradient (stream filament) method and the finite -difference solution of the stream- 
function equation. Each has advantages and limitations. In particular, the finite- 



difference solution of the stream -function equation (e.g. , ref. 1) is limited to strictly 
subsonic flows. The velocity-gradient methods are not limited in this way (e.g. , ref. 2). 
On the other hand, a simple velocity-gradient method is limited to a well-guided channel. 
The purpose of the program described herein is to combine these methods so as to ex- 
tend the range of cases which can be solved. Locally supersonic (transonic) solutions 
can be obtained even with low-solidity blading (channel not well guided). This program 
is called TSONIC (for transonic). 

The TSONIC program is based on the TURBLE program (ref. 3), with the addition 
of subroutines for solving the velocity-gradient equation. Therefore, the input for 
TSONIC is identical to that for TURBLE, but with an additional input item to give a re- 
duced weight flow factor which is needed to obtain a preliminary subsonic solution. 

TSONIC obtains the numerical solution for ideal, transonic, compressible flow for 
an axial, radial, or mixed flow cascade of turbomachine blades. The cascade may be 
circular or straight (infinite) and may be fixed or rotating. To accommodate either 
radial or axial flow with the same coordinate system, the independent variables are 
meridional streamline distance and angle in radians. 

This report includes the FORTRAN IV program TSONIC with a complete description 
of the input required. Numerical examples have been included to illustrate the use of 
the program. Only the parts of the program which differ from the TURBLE program are 
described, although the complete program listing is given. 

This report is organized so that the engineer desiring to use the program needs to 
read only the sections MATHEMATICAL ANALYSIS, NUMERICAL EXAMPLES, and 
DESCRIPTION OF INPUT AND OUTPUT. Information of interest to a programmer is 
contained in the sections DESCRIPTION OF INPUT AND OUTPUT and PROGRAM PRO- 
CEDURE. 

A TSONIC source deck on tape is available from COSMIC (Computer Software Man- 
agement and Information Center) , Computer Center, University of Georgia 30601. The 
program can be ordered by using the number of this report as identification. 


SYMBOLS 


A coefficient in differential eq. (4) 

B coefficient in differential eq. (4) 

b stream- channel thickness normal to meridional streamline, meters 
m meridional streamline distance, meters, see figs. 2 and 3 
R gas constant, J/(kg)(K) 

r radius from axis of rotation to meridional stream- channel mean line, meters 
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S streamline distance, meters 

s angular blade spacing or pitch, rad 

T temperature, K 

t time, sec 

u stream function 

V absolute fluid velocity, meters/sec 

W fluid velocity relative to blade, meters/sec 
w mass flow per blade flowing through stream channel, kg/sec 
z axial coordinate, meters 

a angle between meridional streamline and axis of rotation, rad, see figs. 1 and 3 
/3 angle between relative velocity vector and meridional plane, rad, see fig. 1 
y specific-heat ratio 

77 outer normal to region 

0 relative angular coordinate, rad, see figs. 1 and 2 

2 

X prerotation, (rVo) , meters /sec 

in 

3 

p density, kg/meters 

w rotational speed, rad/sec 

Subscripts: 

cr critical velocity 

in inlet or upstream 

j dummy variable 

le leading edge 

m component in direction of meridional streamline 

out outlet or downstream 

r radial component 

te trailing edge 

z axial component 

0 tangential component 
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Superscripts: 

* absolute stagnation condition 

" relative stagnation condition 


MATHEMATICAL ANALYSIS 

The calculations are performed in two stages. The first stage is to obtain a solu- 
tion based on a reduced weight flow by the finite-difference solution of the stream- 
function equation as described in references 1 and 3. For this first stage of the calcu- 
lations, weight flow must be reduced sufficiently so that the flow is completely subsonic 
throughout the passage. 

The second stage is to obtain a velocity distribution based on the actual weight flow 
by means of a velocity-gradient equation. The velocity-gradient solution requires infor- 
mation obtained in the first stage. There may be locally supersonic flow. 

The velocity across the width of a curved passage will vary. Hence, at the throat 
of a curved passage that is choked, there will be both supersonic and subsonic velocities 
across the passage width. If the weight flow is just slightly less than choking, two solu- 
tions are possible, and both solutions will have both subsonic and supersonic velocities. 
However, the solutions obtained by TSONIC are always the overall "subsonic" solution 
(i.e . , the velocities are always less than those corresponding to choking weight flow). 

The simplifying assumptions used are those in references 1 and 3. These assump- 
tions are 

(1) The flow is steady relative to the blade. 

(2) The fluid is a perfect gas (constant c ) or is incompressible. 

(3) The fluid is nonviscous, and there is no heat transfer (therefore, the flow is 
isentropic). 

(4) The flow is absolutely irrotational. 

(5) The blade -to -blade surface is a surface of revolution. (This does not exclude 
straight infinite cascades.) 

(6) The velocity component normal to the blade-to-blade surface is zero. 

(7) The stagnation temperature is uniform across the inlet. 

(8) The velocity magnitude and direction are uniform across both the upstream and 
downstream boundaries. 

(9) The only forces are those due to momentum and pressure gradient. 

These assumptions are used throughout the program. The following assumption is used 
only in the first stage of the calculation: The relative velocity is subsonic everywhere. 
(This is accomplished by using a reduced weight flow in the first stage of the calculation.) 
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The flow may be axial, radial, or mixed and there may be a variation in the normal 
stream-channel thickness b in the through-flow direction. Since the stream- channel 
thickness can be specified as desired, a loss in relative stagnation pressure can be 
accounted for by reducing b by a percentage equal to the percentage loss in relative 
stagnation pressure. 

The notation for velocity components is shown in figure 1. For generality, the 



Figure L - Cylindrical coordinate system and velocity components. 


meridional streamline distance m is used as an independent coordinate (see fig. 2). 
Thus, m and 6 are the two basic independent variables. A stream channel is defined 
by specifying a meridional streamline radius r and a stream-channel thickness b as 
functions of m alone (see fig. 3). The variables r and b are constant functions of 6 . 
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Figure 2. - Biade-to-blade surface of revolution, 
showing m-0 coordinates. 



1 1 

Figure 3. - Flow in a mixed-flow stream channel. 


For the mathematical formulation of the problem the stream function is used. The 
stream function u is normalized so that u is 0 on the upper surface of the lower 
blade, and 1 on the lower surface of the upper blade. The stream function satisfies the 
following equation (ref. 1). 


:5i: 
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( 1 ) 


1 3^u + 3^u 1 1 3p 3u | pin a _ _1_ 3(bpjj 3u _ 2bpo; g ^ n a 

r 2 59^ 3m^ r 2 p 39 39 L r bp 3mJ 3m w 

The derivatives of the stream function satisfy 

J± = -*£ W fl 
3m w 


( 2 ) 


iH = S££w m (3) 

30 W m 



If the flow is entirely subsonic, equation (1) is elliptic. Boundary conditions for the 
entire boundary ABCDEFGHA of figure 4 will determine a unique solution for u. These 
boundary conditions (ref. 1) are as follows: 
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Boundary segment 

Boundary condition 

AB 

u is 1 less than the value of u on GH at 
the same m coordinate. 

BC 

u = 0 

CD 

u is 1 less than the value of u on EF at 
the same m coordinate. 

DE 

© 

. _ tan/3 out 
out sr out 

EF 

u is 1 greater than the value of u on CD at 
the same m coordinate. 

FG 

u = 1 

GH 

u is 1 greater than the value of u on AB at 
the same m coordinate. 

AH 

© 

_ tan p in 
in sr in 


For the case where there is locally supersonic flow, equation (1) is no longer ellip- 
tic in the entire region, but is hyperbolic in the region of supersonic flow. (This is 
discussed in chapter 14 and appendix A of ref. 4.) With a mixed-type problem like 
this, an analytical solution to equation (1) probably does not exist. This is discussed in 
reference 5 and means simply that there is probably a shock loss. However, the shock 
loss may be so small as to be negligible in a numerical solution. In this case we are 
justified in looking for a numerical solution to equation (1). 

At first one may think that equation (1) could be solved by a finite-difference method 
even when there is locally supersonic flow. There are, however, difficulties with this 
approach. The difficulty has to do with the fact that there are two velocities, one sub- 
sonic and one supersonic, which will give the same value for the weight flow parameter 
pW. If a stream-function solution is obtained, we can calculate the stream-function 
derivatives to obtain values of pW by using equations (2) and (3). However, if there is 
locally supersonic flow, there is no easy way of telling which points should use the sub- 
sonic velocity and which points should use the supersonic velocity. This is further com- 
plicated by the fact that equation (1) is nonlinear and requires iteration to obtain the coef- 
ficients involving the density p. Therefore, in the initial iteration the predicted values 
of pW near the supersonic region will be too large, so that no velocity W can be 
found to correspond to the predicted value of pW. 

Because of the difficulties with the finite -difference method of solution a different 


m 

m 
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technique is used for the case with locally supersonic flow. The method is based on the 
following velocity-gradient equation: 


= AW + B 

00 


( 4 ) 


where 


A = r 2 cos 2 B + sin a tan ^ (1 + cos 2 B) 
dm 2 


(5a) 


is used on blade surface , 


A = sin 2 B 


a 2 
0 u 

dd dm 

d u 

dm 


du 

de 


0 2 u 


' 8u \ 2 0m 2 
0m/ 


a 2 " 

0 U 


d<r 

du 

00 


+ sin a tan B (1 + cos 0) 


(5b) 


is used at interior points , and 


x, „ 0W . 2u>r sin a 

B = r tan B + 

0m cos B 


( 6 ) 


Equations (4) to (6) are derived in appendix A from the force equation. The quan- 
tities in equations (5) and (6) are known if a solution to equation (1) is known. Therefore, 
it is desired to obtain an approximate solution to equation (1) at a reduced weight flow 
such that the streamlines at the reduced weight flow correspond closely to those for the 
actual weight flow. If the flow were incompressible, there would be no change at all in 
the streamlines if w is reduced by the same ratio as the weight flow w. This can be 
seen from equation (1), since the coefficients are constants for incompressible flow and 
equation (1) does not change if co/w is kept constant. None of the quantities in equa- 
tions (5) and (6) would change, except for 0W/0m, which is proportional to the weight 
flow w for incompressible flow. For the compressible case there will be some change 
in streamline shape, and 0W/0m will not be strictly proportional to the weight flow. 
However, for many cases, the error introduced by using quantities in equations (5) 
and (6) from a reduced weight flow solution is not significant. Hence, a solution to 
equation (1) can be obtained for a reduced weight flow (with correspondingly reduced u>), 
such that the values of fi, 0W/0m, and the partial derivatives of u can all be estimated 
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reasonably. More complete detail on obtaining the reduced weight flow solution is given 
in appendix B. Alternate expressions for the parameter A are given since d^0/dm^ 

is known on the blade surface, but the partials of u are known in the interior of the 
region. 

Equation (4) can be solved numerically along vertical line in figure 4 (i.e. , con- 
stant m) by using the approximate values of A and B calculated from the completely 
subsonic solution. To solve equation (4) directly requires a known value of W at some 
initial point, for example, the intersection of the vertical line with the lower boundary of 
the region. However, the condition that determines a unique solution to equation (4) is 

not a value of W at some point, but rather that continuity must be satisfied. That is, 
we require that 


/ 2 

p W cos /3 br d0 = w (7) 

where 9^ is the value of 9 at the lower boundary and 9^ is the value of 9 at the upper 
boundary. One way to satisfy equation (7) is to try some initial value of W, and then 
solve equation (4), and calculate the integral in equation (7). If the value of the integral 
is not equal to w, then other initial values can be tried successively, iterating until the 
correct solution is obtained. The numerical solution of equation (4) is obtained by a 
Runge-Kutta method as described in the description of subroutine VELGRA. 

When equation (4) has been solved, subject to satisfying equation (7), the velocity 
distribution is known along vertical line running from the lower to the upper boundary of 
the region (fig. 4). By doing this for a large number of vertical lines, we obtain the 
velocity distribution for the entire region, including both blade surfaces. 


NUMERICAL EXAMPLES 

To illustrate the use of the program and to show the type of results which can be 
obtained, two numerical examples are given. The first example is an axial flow turbine 
stator, and the other is a radial inflow turbine. 


Axial Stator 

This example is a stator nozzle mean blade section (fig. 5) for a turbine built at 
Lewis (ref. 6). This blade section has also been analyzed by using the subsonic com- 
pressible flow program TURBLE of reference 3. These results are the same as that 
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Figure 5. - Axial stator blade for numerical example. 


reported in reference 7. The design weight flow could not be analyzed with the TURBLE 
program, because the velocities were too close to the sonic velocity. However, satis- 
factory results were obtained with TSONIC by using a value of REDFAC = 0. 8. This 
means that the stream-function solution was obtained first with 80 percent of design 
weight flow, then the velocity-gradient method was used to obtain the velocities at design 

weight flow. 

The input for this case is given in table I. The design weight flow velocities calcu- 
lated by the program are given in table n. These velocities are plotted against blade 
surface length in figure 6. Also shown in figure 6 are experimental data obtained from 
the investigation described in reference 6. The experimental velocities shown are taken 
from figure 13(b) of reference 6. Most of the velocities are in very close agreement. 
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Figure 6. - Blade surface velocities for axial stator mean 
section compared with experimental data. 


Radial Inflow Turbine 

An example turbine rotor profile is shown in figure 7. This is a high-specific- 
speed turbine, with a specific speed of 1.01 (131 (rpm)(ft)^/^/(sec)*/^). First a meri- 
dional plane analysis was made by the quasi -orthogonal method described in reference 11. 
This analysis determined the meridional streamline spacing. The meridional stream- 
line spacing gives the information for array BESP in the input. A preliminary solution 
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is presented here for a blade -to -blade stream channel adjacent to the shroud. The 
outlet flow angle /3^ e is estimated based on the blade trailing- edge angle. The solution 
will help decide whether the estimated value of /3^ e is correct. 

The blade -to -blade channel is shown in figure 8. The input data for the program 


H G 



are given in table HE. The final velocities, based on the full weight flow, are given in 
table IV. These velocities are plotted against the meridional streamline distance m in 
figure 9. 

Figure 9 indicates a high blade loading with moderate suction surface diffusion. 
However, the blade loading probably is higher than would actually be obtained near the 
trailing edge. This indicates that, most likely, less turning would be obtained than that 



Figure 9. - Blade surface velocities for radial turbine 
rotor. 
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TABLE in. - INPUT FOR RADIAL INFLOW TURBINE CASE 
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1 


given as input to the program. The curves for the suction and pressure surface veloc- 
ities would be made to come together at the trailing edge by reducing ft in additional 
computer runs. e 


DESCRIPTION OF INPUT AND OUTPUT 

The computer program requires as input a geometrical description in m-0 coor- 
dinates of the blade surfaces, a description in m-r coordinates of the stream channel 
through the blades, appropriate gas constants, and operating conditions such as inlet 
temperature and density, inlet and outlet flow angles, weight flow, and rotational speed. 
Figures 1 and 2 show the m-0 coordinate system for a typical blade-to-blade surface 
of revolution. Output obtained from the program includes velocity magnitude and direction 
at all interior mesh points in the blade-to-blade passage, blade surface velocities, 
stream-function values throughout the blade-to-blade region of solution, and streamline 
locations . 


1 Jj6 join 15| 16 20j21 25126 30|31 3513* 40|41 50|51 fflfo] 7 nlri an 
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Figure 10. - Input form. 
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Input 


Figure 10 shows the input variables as they are punched on the data cards. There 
are two types of variables, geometric and nongeometric. The geometric input variables 
are shown in figures 11 and 12. All input variables are described in this section. 



Figure 11. - Geometric input variables on blade-to-blade stream surface. Angles BETA!, BETAO, BETH, BETI2, 
BET01, and BET02 must be given as true angle 3 in degrees, not the angle as measured in m-0 plane. 

Either use tan 3 = r d6/dm to obtain 3> or measure the true angle. 



CD-10241 o Spline point 


Figure 12. - Geometric input variables describing stream channel in meridional plane. 


Further explanation of key variables is given in the section Instructions for Preparing 
Input. 

The input variables are as follows: 

GAM specific-heat ratio 
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AR 

TIP 

RHOIP 

WTFL 

OMEGA 

ORF 

BETAI 

BETAO 

CHORDF 

STGRF 

REDFAC 

DENTOL 


MBI 

MBO 

MM 

NBBI 

NBL 

NRSP 


gas constant, J/(kg)(K) 

inlet stagnation temperature, K 

inlet stagnation density, kg/meter 3 

mass flow per blade for stream channel, kg/sec 

rotational speed, o>, rad/sec (Note that w is negative if rotation is in the 
opposite direction of that shown in fig. 1.) 

value of overrelaxation factor to be used in the solution of the inner iteration 
simultaneous equations (If ORF = 0, the program calculates an estimated 
value for the overrelaxation factor. See p. 25 for discussion.) 

inlet flow angle /3j g along BG with respect to m-direction, deg, see fig. 11 

outlet flow angle /3^ e along CF with respect to m-direction, deg, see 
fig. 11 

overall length of blade in m-direction, meters, see fig. 11 

angular d - coordinate for center of trailing-edge circle of blade with respect 
to the center of leading-edge circle of blade, radians, see fig. 11 

factor by which weight flow (WTFL) must be reduced in order to assure sub- 
sonic flow throughout passage (REDFAC is usually between 0. 5 and 0.9.) 

tolerance on density change per iteration for reduced weight flow (DENTOL 
may be left blank, and the value 0.001 will be used. If trouble is exper- 
ienced in obtaining convergence (i.e. , the maximum relative change in 
density (item 14 of the output) does not get small enough), then a larger 
value of DENTOL may be used, or a smaller value of REDFAC may be 
used. The value of 0. 001 for DENTOL is a tight tolerance, 0. 01 would be 
a medium tolerance, and 0. 1 would be a loose tolerance. ) 

number of vertical mesh lines from AH to BG inclusive, see fig. 13 

number of vertical mesh lines from AH to CF inclusive, see fig. 13 

total number of vertical mesh lines in m-direction from AH to DE, maxi- 
mum of 100, see fig. 13 

number of mesh spaces in 0-direction between AB and GH, maximum of 50, 
see fig. 13 

number of blades 

number of spline points for stream- channel radius (RMSP) and thickness 
(BESP) coordinates, maximum of 50, see fig. 12 
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RI1,RI2 
R01,R02 
BETI1, BETI2 


BET01,BET02 


SPLN01,SPLN02 


MSP1,MSP2 


THSP1,THSP2 


MR 


RMSP 

BESP 


leading-edge radii of the two blade surfaces, meters, see fig. 11 

trailing-edge radii of the two blade surfaces, meters, see fig. 11 

angles (with respect to m- direction) at tangent points of leading- 
edge radii with the two blade surfaces, deg, see fig. 11 (These 
must be true angles in degrees . If angles are measured in the 
m-0 plane, (i.e. , d0/dm), BE Til and BETI2 can be obtained 
from the relation tan /3 = r(d0/dm).) 

angles (with respect to m-direction) at tangent points of trailing- 
edge radii with the two blade surfaces, deg, see fig. 11 (These 
must also be true angles in degrees, like BETI1 and BETI2.) 

number of blade spline points given for each surface as input, max- 
imum of 50 (These include the first and last points (dummies) 
that are tangent to the leading- and trailing-edge radii (fig. 11).) 

arrays of m-coordinates of spline points on the two blade surfaces, 
measured from the blade leading edge, meters, see fig. 11 (The 
first and last points in each of these arrays can be blank or have 
a dummy value, since these points are calculated by the program. 
If blanks are used, and the last point is on a new card, a blank 
card must be used.) 

arrays of 0-coordinates of spline points corresponding to MSP1 
and MSP2, radians, see fig. 11 (Dummy values are also used 
here in positions corresponding to those in MSP1 and MSP2. ) 

array of m-coordinates of spline points for stream- channel radii 
and stream- channel thickness, meters, see fig. 12 (MR is 
measured from the leading edge of the blade. These coordinates 
should cover the entire distance from AH to DE, and may extend 
beyond these bounds. The total number of points is NRSP.) 

array of r-coordinates of spline points for the stream -channel radii, 
corresponding to the MR array, meters, see fig. 12 

array of stream -channel normal thicknesses corresponding to the 
MR and RMSP arrays, meters, see fig. 12 


The remaining variables, starting with BLDAT, are used to indicate what output is 
desired. A value of 0 for any of these variables will cause the output associated with 
that variable to be omitted. A value of 1 will cause the corresponding output to be 
printed for the final iteration only; 2, for the first and final iterations; and 3, for all 
iterations. Care should be used not to call for more output than is really useful. The 
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following list gives the output associated with each of these variables: 

BLDAT all geometrical information which does not change from iteration to iteration 
(i.e. , coordinates and first and second derivatives of all blade surface 
spline points; blade coordinates, blade slopes, and blade curvatures where 
vertical mesh lines meet each blade surface; radii and stream -channel 
thicknesses corresponding to each vertical mesh line; m- coordinate, stream- 
channel radius and thickness, and blade surface angles and slopes where 
horizontal mesh lines intersect each blade; and ITV and IV arrays, internal 
variables describing the location of the blade surfaces with respect to the 
finite difference grid) 

AANDK coefficient array, constant vector, and indexes of all adjacent points for each 
point in finite -difference mesh (This information is needed for debugging 
the program only. ) 

ERSOR maximum change in stream function at any point for each iteration of SOR 
equation (eq. (A8), ref. 1) 

STRFN value of stream function at each unknown mesh point in region 

SLCRD streamline 0- coordinates at each vertical mesh line, and streamline plot 

INTVL velocity and flow angle at each interior mesh point for both reduced and actual 
weight flow 

SURVL m-coordinate, surface velocity, flow angle, distance along surface, and 

W/W cr based on meridional velocity components where each vertical mesh 
line meets each blade surface; m-coordinate, surface velocity, flow angle, 
distance along surface, and W/W cr based on tangential velocity components 
where each horizontal mesh line meets each blade surface; plot of blade 
surface velocities against meridional streamline distance, meters 


Instructions for Preparing Input 

It is very unusual to have no errors of input the first time TSONIC is run. There- 
fore, it is recommended that the first attempt should allow only 1 minute of execution 
time and that BLDAT should be equal to 1. The resulting output should be checked care- 
fully. Of particular interest are the second derivatives at input spline points. Any 
errors in blade geometry input will usually result in wild values for some of these sec- 
ond derivatives. All other preliminary output should be checked to see that it is reason- 
able. 
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Units of measurement. - The International System of Units (ref. 8) is used through- 
out this report. However, the program does not use any constants which depend on the 
system of units being used. Therefore, any consistent set of units may be used in pre- 
paring input for the program. For example, if force, length, temperature, and time 
are chosen independently, mass units are obtained from Force =Mass x Acceleration. 

The gas constant R must then have the units of force times length divided by mass 
times temperature (energy per unit mass per degree temperature). Density is mass per 
unit volume, and weight flow is mass per unit time. Output than gives velocity in the 
chosen units of length per unit time. Since any consistent set of units can be employed, 
the output is not labeled with any units. 

Blade and stream- channel geometry . - The upper and lower surfaces of the blade 
are each defined by specifying three things: leading- and trailing-edge radii, angles at 
which these radii are tangent to the blade surfaces, and m- and B - coordinates of sev- 
eral points along each surface . These angles and coordinates are used to define a cubic 
spline curve fit (ref. 9) to the surface. The standard sign convention is used for angles, 
as indicated in figure 11. 

A cubic spline curve is a piecewise cubic polynomial which expresses mathematically 
the shape taken by an idealized spline passing through the given points. Reference 9 
describes a method for determining the equation of the spline curve. Using this method, 
few points are required to specify most blade shapes accurately, usually no more than 
five or six in addition to the two end points. As a guide, enough points should be spec- 
ified so that a physical spline passing through these points would accurately follow the 
blade shape. This means that the spline points should be closer where there is large 
curvature and farther apart where there is small curvature. As a check, the program 
should be run for 1 minute of execution time with BLDAT = 1 for any new geometry. 

Check the second derivatives at the spline points to see that they are reasonable. Also 

check blade curvatures at vertical mesh lines. 

The coordinates for either surface of the blade are given with respect to the leading 
edge, with the leading edge of the blade being defined as the furthest point upstream (see 

fig- 11)- 

The mean stream surface of revolution (as seen in the meridional plane, fig. 12) and 
the stream- channel thickness are also fitted with cubic spline curves. The m- 
coordinates for the mean stream surface are independent of the m-coordinates for the 

blade surface. 

Loss of relative stagnation pressure . - If desired, a simplified correction for losses 
can be made by assuming a loss in relative stagnation pressure. This type of loss can 
be accounted for by reducing each value in the BESP array by a percentage equal to the 
assumed percentage loss in relative stagnation pressure at that point. 

Inlet and outlet flow angles . - The values of 0 le and 0 te are given as average 
values on BG and CF, respectively. If the flow is axial, these flow angles are constant 
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upstream or downstream of the blades. If flow is radial or mixed, and these angles are 
not known on BG and CF, /3j e and /3 tfi must be calculated by equation (B14) of refer- 
ence 1 or (B15) of reference 7. 

Defining the mesh . - A finite-difference mesh is used for the solution of the basic 
differential equation. A typical mesh pattern is shown in figure 13. The mesh spacing 
and the extent of the upstream and downstream regions are determined by the values of 
MBI, MBO, and MM of the input. The mesh spacing must be chosen so that there are 
not more than 2500 unknown mesh points. 



Values of MBI, MBO, and MM should be determined so that the mesh which results 
has blocks which are approximately square. To achieve this, a value for NBBI is first 
chosen arbitrarily (15 to 20 is typical). NBBI is the number of mesh spaces spanning 
the blade pitch s, where s = 2 tt/NBL. Dividing s by NBBI gives the mesh spacing HT 



1*1 


= 8-5 
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in the 0-direction in radians. Multiplying HT by an average radius (RMSP) of the stream 
channel gives an average value for the actual mesh spacing in the 0-direction. The 
value of CHORD should then be used with this tangential mesh spacing to calculate the 
approximate number of mesh spaces along the blade in the m-direction. This will give 
MBO once MBI is chosen. Generally, MBI is given a value of 10. MM, likewise, is 
usually given a value 10 more than MBO. 

Over relaxation factor. - ORF is the over relaxation factor used in each inner itera- 
tion in the solution of the simultaneous finite -difference equations (see ref. 2, p. 102). 
ORF may be set to 0, or some value between 1 and 2. ORF should be 0 for the initial 
run of a given blade geometry and mesh spacing (MBI, NBBI, etc.). In this case the 
program uses extra time and calculates an optimum value for ORF . It does this by 
means of an iterative process, and on each iteration the current estimate of the optimum 
value for ORF is printed. The final estimate is the one used by the program for ORF . 

If the user does not change the mesh indexes MBI, MBO, MM, and NBBI between runs, 
even though blade geometry or other input does change, he may use this final estimate 
of ORF in the input, saving the time used in its computation. In all cases, if ORF is not 
0, it should have a value greater than 1 and less than 2. 

Actually, the value of ORF is not as critical as the user might think. It gets more 
critical as the optimum value gets close to 2. For any run of a given set of data, only 
small changes will occur in the rate of convergence in SOR as long as the difference 
2.0 - ORF is within 10 percent of its optimum value. 

Format for input data . - All the numbers on the card beginning with MBI and on the 
card beginning with BLDAT are integers (no decimal point) in a five-column field (see 
fig. 10). These must all be right adjusted. The input variables on all other data cards 
are real numbers (punch decimal point) in a 10- column field. 

Incompressible flow. - While the program is written for compressible flow, it can 
be easily used for incompressible flow. To do so specify GAM = 1.5, AR = 1000, 

TIP = 10 6 , and REDFAC = 1 as input. This results in a single outer iteration of the pro- 
gram to obtain the stream -function solution. And, of course, the velocity-gradient solu- 
tion will yield nothing new. 

Straight infinite cascade . - The program is as easily applied to straight infinite cas- 
cades as to circular cascades. Since the radius and number of blades (NBL) for such a 
cascade would actually be infinite, an artificial convention must be adopted. The user 
should pick a value for NBL, for instance 20 or 30. Then, since the blade pitch sr is 
known, an artificial radius can be computed from 

r _ NBL(sr) 

277 

This radius should then be used to compute the 0-coordinates required as input (THSP1, 
THSP2 , and STGRF). 
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Axial flow with constant stream -channel thickness . - For a two-dimensional cascade 
with constant stream -channel thickness, constant values should be given for the MR, 
RMSP, and BESP arrays. Only two points are required for each of these arrays in this 
case. The two values of MR should be chosen so that they are further upstream and 
downstream than the boundaries AH and DE. The two values of RMSP and BESP should 
equal the constants r and b. 


Output 

Sample output is given in table V for the axial flow stator example of reference 1. 
The blade shape is shown in figure 13. Since the complete output would be lengthy, only 
the first few lines of each section of output are reproduced here. Most of the output is 
optional, and is controlled by the final input card, as already described. In some in- 
stances output lables are simply internal variable names. 

Each section of the sample output in table V has been numbered to correspond to the 
following description: 

(1) The first output is a listing of the input data. All items are labeled as on the 
input from (fig. 10). 

(2) This is the output corresponding to BLDAT (see the list of input variables and the 
descriptions of internal variables for the subroutines of the program). 

(3) The relative free-stream velocity W; the relative critical velocity W cr ; and 
the maximum value of the mass flow parameter pW (corresponding to W = w c p are 
given at the leading edge of the blade (BG) and the trailing edge of the blade (CF). These 
are all for the full weight flow. The inlet (outlet) free-stream flow angle /3^ n (/3 qu j.) at 
boundary AH (DE) is given. These angles are based on the reduced weight flow and the 
input angles BETAI (/3j g ) and BETAO (j3 te ). The reason for this is discussed in appen- 
dix B. 

(4) These are calculated program constants, including the pitch from blade to blade, 
the mesh spacing in all solution regions, the minimum and maximum values of IT in the 
solution region (ITMIN and ITMAX), and the value of the prerotation X (eq. (B8), ref. 1) 
for both full and reduced weight flow. 

(5) This is the number of mesh points in the entire solution region at which the 
stream function is unknown. 

(6) This is the boundary value BV of the stream function on each of the blade sur- 
faces . 

(7) This is the output corresponding to AANDK. 

(8) If the program calculates on optimum overrelaxation factor (i.e. , ORF = 0 in the 
input), the successive estimates to the optimum value of the ORF are printed. The last 
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TABLE V. - SAMPLE OUTPUT FOR AXIAL FLOW TURBINE STATOR CASE 
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TABLE V. - Continued. SAMPLE OUTPUT FOR AXIAL FLOW TURBINE STATOR CASE 
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TABLE V. - Continued. SAMPLE OUTPUT FOR AXIAL FLOW TURBINE STATOR CASE 
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TABLE V. - Concluded. SAMPLE OUTPUT FOR AXIAL FLOW TURBINE STATOR CASE 


in 

in CD 
o <r 
• • 
in r>» 
r» vO 


m r» 
m co 
o O' 

■ i 

m r- 

r» *o 


in — < 

—* CD 
• • 
in r~ 
r- o 




36 


bLADE SURFACE VELOCITIES 
(bASED ON FULL WEIGHT FLOW) 
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printed value of the estimate optimum ORF is the value of the over relaxation factor (ORF) 
used by the program. 

(9) This is the output corresponding to ERSOR. 

(10) This is the output corresponding to STRFN. 

(11) This is the total execution time after obtaining the stream function solution for 
each outer iteration. 

(12) This is the output corresponding to SLCRD. 

(13) This is the output corresponding to INTVL for the reduced weight flow. 

(14) This gives the maximum relative change in the density, for each outer iteration. 

(15) This is the output corresponding to SURVL for the reduced weight flow. 

(16) This is the total execution time after all calculations are completed for an outer 
iteration with reduced weight flow. 

Most of the previously described output has been for the reduced weight flow. The 
following output is for the actual weight flow. 

(17) This is the output corresponding to INTVL for the full weight flow. 

(18) This is the output corresponding to SURVL for the full weight flow. 

ERROR CONDITIONS 

1. SPLINT USED FOR EXTRAPOLATION. EXTRAPOLATED VALUE = X.XXX 

SPLINT is normally used for interpolation, but may be used for extrapolation in 
some cases. When this occurs the above message is printed, as well as the input and 
output of SPLINT. Calculations proceed normally after this printout. 

2. BLCD CALL NO. XX 
M COORDINATE IS NOT WITHIN BLADE 

This message is printed by subroutine BLCD if the M- coordinate given this sub- 
routine as input is not within the bounds of the blade surface for which BLCD is called. 

The value of m and the blade surface number are also printed when this happens. This 
may be caused by an error in the integer input items for the program. 

The location of the error in the main program is given by means of BLCD CALL 
NO. XX, which corresponds to locations noted by comment cards at each MHORIZ, 

ROOT, and BLCD call in the program. 
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3. ROOT CALL NO. XX 


ROOT HAS FAILED TO OBTAIN A VALID ROOT 

This message is printed by subroutine ROOT if a root cannot be located. The input 
to ROOT is also printed. The user should thoroughly check the input to the main pro- 
gram. 

The location of the error in the main program is given by means of ROOT CALL 
NO. XX, which corresponds to locations noted by comment cards at each MHORIZ and 
ROOT call in the program. 


4. DENSTY CALL NO. XX 
NER(l) = XX 

RHO*W IS X.XXXX TIMES THE MAXIMUM VALUE FOR RHO*W 

This message is printed if the value of pW at some mesh point is so large that there 
is no solution for the value of p and W. This indicates a locally supersonic condition, 
which can be eliminated by decreasing REDFAC in the input. 

If RHO*W is too large, TSONIC still attempts to calculate a solution. This often 
permits an approximate solution to be obtained which is valid at all the subsonic points 
in the region. In other cases the value of W is reduced at some of the points in ques- 
tion during later iterations, resulting in a valid final solution for these points. The pro- 
gram counts the number of times supersonic flow has been located at any point during a 
given run (NER(l)). When NER(l) = 50, the program is stopped. 

The location of the error in the main program is given by means of DENSTY CALL 
NO. XX, which corresponds to locations noted by comment cards at each DENSTY call 
in the program . 


5. MM, NBBI, NRSP, OR SOME SPLNO IS TOO LARGE 

If this message is printed, reduce the appropriate inputs to their allotted maximum 
values . 

6. INPUT WEIGHT FLOW (WTFL) IS TOO LARGE AT BLADE LEADING EDGE 

This message is printed if WTFL is greater than the choking mass flow for the ver- 
tical line BG, and the program is stopped. If this happens, there is probably an error in 
the input. The following items should be checked carefully: RHOIP, WTFL, BETAI, 
NBL, RMSP, and BESP. 
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7. REDUCED WEIGHT FLOW IS STILL TOO LARGE 

This message is printed if difficulty is encountered in calculating /3. or /3 for 
the reduced weight flow. If this happens, REDFAC should be reduced. “ ° Ut 

8. ONE OF THE MH ARRAYS IS TOO LARGE 

This message is printed if there are more than 100 intersections of horizontal mesh 
lines with any blade surface. In this case NBBI should be reduced. 

9. THE NUMBER OF INTERIOR MESH POINTS EXCEEDS 2500 

This message is printed if there are more than the allowable number of finite- 
difference grid points. Either MM or NBBI must be reduced. 

10. SEARCH CANNOT FIND M IN THE MH ARRAY 

If this message is printed, the value of m and the blade surface number are also 
printed. The user should thoroughly check the input to the main program. 

11. A VELOCITY- GRADIENT SOLUTION CANNOT BE 
OBTAINED FOR VERTICAL LINE IM = XX 

This message is printed if difficulty is encountered in solving the velocity- gradient 
equation for some vertical line. 

12. A VELOCITY-GRADIENT SOLUTION COULD NOT BE OBTAINED IN 
50 ITERATIONS FOR VERTICAL LINE IM = XX 

This message is printed after 50 attempts to find a velocity-gradient equation which 
results either in the specified weight flow (WTFL) or in a choked flow. 

13. WTFL EXCEEDS CHOKING WEIGHT FLOW FOR IM = XX 
CHOKING WEIGHT FLOW = XXXXX FOR IM = XX 

This message is printed if the vertical line IM will not pass the specified weight 
flow (WTFL). WTFL should be reduced in this case. 
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PROGRAM PROCEDURE 


The first part of the program is very similar to TURBLE (ref. 3). The program 
description for TURBLE is given in reference 1. The main difference in this part of 
the TSONIC is the calculation of coefficients A and B of equations (4) to (6) by sub- 
routine TANG. Also, PRECAL has been considerably changed to calculate certain con 
stants at reduced weight flow. In addition, a new segment was added to solve the 
velocity- gradient equation. The main subroutine of this new segment is TVELCY. 
PRECAL, TANG, TVELCY and the subroutines in the new segment of the program are 
described later in this section. All the subroutines and their relation are shown in 
figure 14. 



Figure 14. - Calling relation of subroutines. 


The dictionary defines all new variables. Any variables not defined here are de- 
fined in reference 1. 

The program can handle up to 2500 mesh points on the IBM 7094-2/7044 direct 
coupled system with a 32 768- word core. To be able to handle 2500 mesh points, an 
overlay arrangement is used, as shown in figure 15. All subroutines not shown are in 
the main link. If there is a storage problem on the user's computer, the maximum 
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Figure 15. - Arrangement for overlay, shewing octal storage requirements. 


number of mesh points should be reduced. The following program changes are required 
to change the maximum number of mesh points: 

(1) Change the dimension of A, U, K, and RHO in the COMMON/AUKRHO/ state- 
ment. This statement occurs in most subroutines. 

(2) In subroutine INPUT, change the number of values of U, K, and RHO to be 
initialized (the bound on the DO loop near statement 60). 

(3) In subroutine PRECAL, change statement 340 and format statement 1150 to re- 
flect the maximum allowable number of mesh points. Statement 340 will cause the pro- 
gram to stop if there are too many mesh points . 

(4) Change the dimensions of W, BETA, DUDT, DUDTT, AAP, and BBP in SLAX, 
SLAV, TANG, VELOCY, VEL, TVELCY, and VELGRA. 

The conventions used in the program and most of the labeled COMMON blocks are 
described in reference 1. The following are new COMMON blocks in TSONIC: 
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(1) /SURVEL/ is used for arrays of surface velocities and plotting arrays. 

(2) /D2TDM2/ is used for an array of second derivatives d 0/dm along the blade 

surfaces 

(3) /WWCRM/ is used for an array of critical velocity ratios and an array of labels 
to indicate choking. 


Subroutine PRECAL 

Calculation of X, W le , andW te . - The pre rotation X and the average velocity W le 

at the blade leading edge are calculated for the full weight flow. They are calculated 
iteratively by equations (B7) to (B9) of reference 1. H the input weight flow (WTFL) is 
too large, a solution cannot be obtained. In this case, error message (6) is printe , an 
the program is stopped. Otherwise, is calculated using by satisfying one 

dimensional continuity at line CF. 

Calculation of W pr and maximum values of mass flow parameter pW. - Ca cu a ion 
of all these quantities at the leading and trailing edges is done using equations (BIO) to 
(B12) of reference 1. 

Cal culation of w. co, and X for reduced weight flow. - First, the values of w, », 
and X for the full weight flow are stored. Then values of w and co for the reduced 
weight flow are calculated by multiplying by REDFAC. The value of X for the reduced 
weight flow is then calculated by the same procedure as for the full weight flow. 

Calculate 0 in and (3 QuV - These quantities are necessary as boundary conditions. 

The input, however, gives \ e (BETAI) and ^ (BET AO). The desired .values of ^ 

and B are calculated from the input values of 0 le and 0 te , respectively, by using 

equation (B14) of reference 1. The reduced weight flow is used in these calculations for 
reasons explained in appendix B. 

The remaining calculations in PRECAL are the same as for TURBLE, and are 
described in reference 1. 


Subroutine TANG 

Most of this subroutine is the same as in TURBLE, and is described in reference 1. 
After these calculations are complete for a given horizontal mesh line segment, there is 
a check to see if there is convergence of the outer iteration (by checking the value of 
TF.ND) If convergence has been achieved, the coefficients A and B of equation 
calculated. This requires first the calculation of d 2 n/dm 39 and 3W/0m. The calcu- 
lation of 0 2 u/0m 00 is done by using subroutine SPLINE to calculate the partial denva- 
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tive with respect to m of 0u/00. The values of 0u/00 have been previously calculated 
in SLAV. The values of 0W/0m are also calculated by subroutine SPLINE. The values 
of A and B are now calculated at each point by using equations (5) and (6). These 
values are stored in the AAP and BBP arrays. 

Subroutine TVELCY 

TVELCY calls VELGRA and VELGRB to solve the velocity-gradient equation along 
each vertical line. After this is done, the blade surface velocity is printed at each ver- 
tical mesh line. Then these velocities are plotted. 


Subroutine VELGRA 


VELGRA has a second entry point. The main entry point is for vertical mesh lines 
upstream or downstream of the blade. The second entry point is called VELGRB and is 
used for vertical mesh lines between the blades. Most calculations are the same for 
either entry point. The variable AORB is used as a switch to indicate the entry point 
and is used where there is some difference in the calculation, such as at the surface of 


After calculating some constants the values of A and B of equation (4) are placed 

in the A 2 and B2 arrays from the AAP and BBP arrays, respectively, for the interior 

points. The values of A and B on the blade surface are calculated by equations (5a) 
and (6). 

After all the values of A and B are calculated, the velocity-gradient equation (4) 
is solved. An initial estimate of W on the lower boundary is available from the reduced 
weight flow solution. The initial velocity is obtained by dividing the reduced weight flow 
value by REDFAC. A numerical solution to equation (4) is calculated by a Runge-Kutta 
method as follows (ref. 10, p. 233). If W is known at the j th point, W. , is calcu- 
lated by the following algorithm. Let J+1 


w ”:='V<A itl w j * +1+ B j+I )( Vl 



( 8 ) 


Then 


W? i +W**. 
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After the second boundary is reached, the solution is checked by calculating the weight 
flow. The weight flow can be calculated by 


w 


est 



pW cos /3 br d0 


( 9 ) 


using the values of W just calculated. The density p for each W is calculated, and 
the value of /3 from the reduced weight flow solution is used. The stream -channel 
thickness b and the radius are constants which are available in the BE and RM arrays, 
respectively. The integral in equation (9) is calculated numerically by subroutine 
INTGRL. After w egt (WTFLES in program) is calculated, subroutine CONTIN is called. 
CONTIN will give a new initial value for W. The entire procedure is then repeated, with 
CONTIN giving a new initial value for W each time, until one of four occurrences: 

(1) WTFLES = w±w/10 5 . 

(2) A maximum value of WTFLES <C w is found. This indicates choking. In this 
case, error message 13 is printed. 

(3) In the calculations W becomes so large that no corresponding density can be 
found. In this case, error message 11 is printed, and the program goes to the next 
vertical line. 

(4) None of the above conditions are met for 50 iterations. In this case, error 
message 12 is printed, and the program goes to the next vertical line. 

After the calculations are completed, the interior velocities are printed, and surface 
velocities are stored to be printed later. 

Subroutine BLOCKD 

This initializes the LABEL array with blanks. 

Main Dictionary 

Most of the FORTRAN variables are the same as those given in reference 1 and will 
not be repeated here. The following list defines all new variables in the previously dis- 
cussed subroutines . 

A2 array of values of coefficient A (eq. (4)) along a vertical mesh line 

AAP array of values of coefficient A (eq. (4)) at all interior mesh points 

ACTLAM value of X based on input value of weight flow w (WTFL) (LAMBDA is 

based on the reduced weight flow for the reduced weight flow calculations.) 
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ACTOMG 

ACTWT 

AORB 

B2 

BBP 

CBETA 

CHOKED 

D2TDM2 

DDT 

DELMAX 

DENTOL 

DUDMM 

DUDTM 

DUDTT 

DWDM 

12 

IND 

LABEL 

NERT 

REDFAC 

RWCB 

SBETA 

SBETA1 

SBETAN 

THETA 


input value of OMEGA (OMEGA is reduced for the reduced weight flow 
calculations . ) 

input value of WTFL (WTFL is reduced for the reduced weight flow cal- 
culations . ) 


switch in VELGRA to indicate entry point (AORB = 1 for main entry point 
and AORB = 2 for entry point VELGRB. ) 

array of values of coefficient B (eq. (4)) along a vertical mesh line 

array of values of coefficient B (eq. (4)) at all interior mesh points 

array of values of cos /3 along a vertical mesh line 

variable storing the word "choked" in Hollerith 

array of values of d?9/dm^ for each blade surface at vertical mesh lines 

array of value of 3u / 59 along a horizontal line 

maximum permitted change of estimated initial value of W 
(DELMAX = W cr /10) 


see input 

array of values of 3^u/ 3m^ along a horizontal mesh line 
array of values of 5 u/ 30 3m along a horizontal line 
array of values of 3^u/ 30^ at all interior mesh points 
array of values of 3W/ 3m along a horizontal mesh line 
temporary index variable in TVELCY 


integer variable controlling logical sequence in CONTIN 

array of labels for A format output to indicate a particular blade surface 
velocity was based on choked weight flow 

temporary storage of a value of NER(l) 

see input 

array of values of pW cos /3 along a vertical mesh line 
sin p 

sin 0 on the upper blade surface at a given vertical line 
sin j3 on the lower blade surface at a given vertical line 
array of 0- coordinates along a vertical mesh line 


TOLERC velocity tolerance for convergence in calculating choking weight flow 
(TOLERC = W/100) 

TORSAL 2o>r sin a 

VT temporary velocity 

WAS W* +1 , eq. (8) 

WASS W**p eq. (8) 

WGRAD array of velocities W calculated by eq. (8) 

■^YXP array of velocities W along a horizontal mesh line in TANG 

WTFLES w e st ca l cu ^ ate ^ by ec l- W 


Program Listing For Subroutines Using Main Dictionary 


COK.^ON $RW f I TtR f I END * LER ( 2 ) » NER ( 2 ) 

COMMON / AUKRHO/ A(2500,4),U(2500),K(2500),RHG(2500) 

COMMON / I NP/GAM , AR »T I P *RKH P »WTFL » OMEGA. ORF , BE TA I , BETAO, REDF AC , 

1 OENTOL.MGI , MBO , MM , NBB 1 , NBL , NRSP . MR ( 50 ) , RMSP ( 50 ) . BFSP ( 50 ) , 

2 BLOAT , AANDK , ERSOR » STRFN, SLCRD, INTVL .SURVL 

COMMON /CALCON/ACTWT » ACTOMGt ACTLAM »MB IM1 » MBI PI »MBOM l »MB0P1 » MMH l » 


/IALUJIM / I W I f i unvi I I ' T _ T T n n lAMDPiA 

HMl, HT.DTLR.OMLR, PITCH, CP, EXPON.ThW.CPTP.TGROG.TBI.TBO. LAMBDA. 

TWL* ITMIN.1T MAX , N I P » I MS ( 2 ) , 3V ( 2 ) » MV ( 100 ) , I V ( 10 1 ) , I TV ( 100 , 2 » 

TV( 100, 2) ,UT0MV( 100, 2 ), BETAVl 100,2).MH(100,2) ,DTDMH I 100, 2 ) , 

BET AM ( 100, 2 ) , RMH( 1C0,2),BEH( 1 00 , 2 ) , RM ( 100 ) , OE ( 1 00 ) , DBDM { 1 00 ) , 

COMMON ^ /GEOMIN/ CHORD l 2 ),STGR12),MLE(2) , THLE ( 2 ) , RM I (2 ) ,RMO ( 2 J , 

1 8 1 ( 2 ) » RO ( 2 ) » BET I ( 2 I »BCT0 ( 2 ) » NSP I ( 2 ) » MSP ( 50, 2 ) » THSP ( 50 » 2 ) 

COMMON /RHOS/RHOHB ( 100 ,2 ) » RHGVBl 100,2) 

COMMON / BLCUCM/ EM ( 50 , 2 ) , I NI T ( 2 ) vironuinm 

COMMON /SURVEL/ WTB( ICO, 2) ,WMB( 100,2), XDOhN (A00) .YACROS(AOO) 

COMMON /D2TDM2/ D2TDM 2 ( 100 , 2 ) ,. T cun cnuc cmCT 

INTEGER BLOAT, AANDK.ER SOR, STRFN, SLCRD, SURVL, AATEMP.SURF.FI RST, 

^ PE AU K.KAK, LAMBDA, LMAX.MP.MLE, MR, MSL.MSP.MV.MVIMl 


CALL TIMElCTll 
10 I END = -1 
ITER = 0 
INIT(l) = 0 
IN I T ( 2 ) = 0 
CALL INPUT 
CALL PKECAL 
30 CALL COEF 
CALL SOR 
CALL T I ME 1 ( T 2 ) 

T I ME= ( T2-T 1 ) / 3600* 
VnRI TE(6, 1000) TIME 
CALL SLAX 
CALL TANG 



u u o o o o u 


CALL VELOCY 

CALL TIMEKT21 

T I ME = (T2-T1 )/3600. 

WRI TEI6, 1000 ) TIME 
I F ( NER ( 2 ) . GT . 0 ) GO TO 10 
IF ( I END . L E . 0 ) GO TO 30 
CALL TVELCY 
GO TO 10 

1000 FORMAT (8HLTIME = ,F7.4,5H MIN.) 
END 


SUBROUTINE INPUT 

INPUT READS ANO PRINTS ALL INPUT DATA CARDS AND CALCULATES HORIZONTAL 
SPACING (MV ARRAY) 

COMMON SR W, ITER, I END , LER (? ) , NER ( 2 ) 

COMMON / AUKKHO/ A ( 250 0 , A J , U { 2500 ) , K < 2500 ) ,RHO ( 2500 ) 

COMMON /INP/GAM,AR,TIP,RhOIP,WTFL,OMEGA,ORF,BtTAI,BETAO»REDFAC, 

1 DENTOL ,MB I , MBO , MM , NBB I , NBL , NR SP , MR ( 50 ) , RMSP ( 50 ) ,BESP( 50) , 

2 BLOAT , AANDK.ERSOR, STRFN.SLCRD, INTVL.SURVL 

COMMON /CALCQN/ACTWT , ACT0MG,ACTLAM,MBIMI,MCIP1 ,MB0M1 ,MB0P1 , MMM 1 , 

1 HMl , HT , OTLR , DMLR, P ITCH, C P , EXPON , TWW , C PT IP , TGROG , TB I ,TBO, LAMBDA, 

2 TWL, ITMIN, IT MAX, NIP, IMS(2 ),BV(2) , MV (100), IV1101) , ITV( 100,2) , 

3 TV ( 100,2) ,DTDMV(100,2 ),BETAV( 100 , 2 ) , MH( 100, 2 ) , DTDMH ( 100,2) , 

A BET AH ( 100,2) ,RMH( 100,2) ,BEH( 100,2),RM( 100 ) , BE ( 1 00 ) , DBDM (100) , 

5 SAL( 100} , AAAI 100) 

COMMON /GEOMIN/ CHOR D ( 2 ) ,STGR ( 2 ) , MLE { 2 ) , THLF (2 ) , RM I ( 2 ) , RMO ( 2 ) , 

1 RI(2) ,R0(2) , BET I( 2 ) ,BETO( 2) ,NSP H 2 ) , MSP I 50, 2 ) , THSP150.2) 

COMMON /RHOS/RHOHB(100,2 ),RHOVB( 100,2) 

INTEGER BLOAT, AANDK, ERSOR, STRFN»SLCRD,SURVL»AA TEMP, SURF, FIRST, 

1 UPPER, SI, ST, SRW 

REAL K,KAK,LAM8DA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1 

REAC AND PRINT ALL INPUT OATA 

WRITEI6.1000) 

READ (5,1100) 

WR I TE ( 6 , 1 100 ) 

WR I TEI 6 , l 1 10 ) 

READ (5,1030) GAM , AR , T I P ,RHG I P , WTFL , BL ANK, OMEGA, ORF 
WRI TEI 6, 1040 ) GAM, AR , T I P ,RHO I P , WTFL , BLANK , OMEGA , ORF 
WRITE(6, 1120) 

READ ( 5 , 1030JBETAI ,BETAO,CHORD( 1 ) ,STGR( 1 ) 

UR I TE( 6, 1040 (BETA I ,BET AO, CHORD ( 1 ) ,STGR( 1 ) 

WRITE (6,1125) 

READ (5,1030) REDFAC , DENTOL 
IF ( DENTOL . LE .0 . ) DENTOL = .001 
WRITE (6,1040) REDFAC .DENTOL 
WRITE16, 1130) 

READ (5,1010) MBI.MBO, BLANK, BLANK, MM, NBBI, NBL, NRSP 
WRITE(6, 1010) MB I , MB 0 , BL ANK , BLANK , MM , NBB I , NBL , NRSP 
CO 10 J=1 , 2 

IF (J.EQ.l) WRITEI6, 1140) 
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non non 


IF (J.EQ.2) WRITE (6, 1150) 

WRITEI6, 1180) J.JtJ.J.J 

READ ( 5 , 1 0 TO ) K I ( J ) , RC ( J ) , BE T II J ) , BE TO ( J ) , SPLNO 
WRI IE ( 6 j 1 0 AO ) RI ( J ) y R G ( J ) » BE T I (J)»BGT0(.J) » SPLNO 
KSP I ( J ) = SPLNO 
NSP = NSP I ( J ) 

WRITE I6»1190) J 

READ (5,1030) (MSPII.J), 1=1, NSP) 

WR I TEC 6 , 1040 ) (MSP ( I , J ) , 1=1, NSP ) 

WRITE! 6, 1200) J 

READ (5,1030) ( THS P ( I , J ) , I = l , NSP ) 

10 WRI IE(6, 10'tO) (THSP( I ,J) ,1 = 1, NSP) 

WRITE(6, 1210) 

RE AO (5,1030) IMR( I) , 1 = 1 ,NRSP) 

WR I TE ( 6 , 1 040 ) ( MR ( I ) , 1=1 ,NRSP) 

WR I TE ( 6 , l 220 ) 

READ (5,1030) ( RMSP( I ) , I =1 , NRSP ) 

WR I TE ( 6 , 1 040 ) (RMSPl I ), 1=1, NRSP) 

WRITE(6,1230) 

READ (5,1030) ( B E S P ( I ), 1=1, NRSP) 

WRITE(6,1040) ( BESPl I ) , I =1 ,NRSP ) 

REAl) E C 5 1 1010 1 BLDAT, A ANDK, ER SOR , STRFN, SLCRD, I NTVL . SURVL 
WRITE16, 1020) BLOAT, AANOK.ERSOR, STRFN, SLCRD, INTVL, SURVL 

IF (MM.LE.100. 4ND.NBB I .L t. 50 . AND. NRSP .LE . 50. AND . NSP I l l ) . LE . 50 
1 .AND.NSPI (2 ) .LE.50) GO TO 20 
WRITE (6,1250) 

STOP 

CALCULATE MV ARRAY 

20 HMl = CHORD! 1 ) /FLOAT ( MB0-M9 I ) 

CO 30 I M= 1 * MM 

30 MV(IM) = FLUAT t IM-MBI ) *HM1 
MV(MBO) = CHOROIl) 

CALCULATE MISCELLANEOUS CONSTANTS 

NER l 1) =0 
NER l 2 ) =0 

PITCH = 2 . *3. 1 415927/ FLOAT ( NBL ) 

HT= PI TCH/FLCaT(NBBI ) 

CTLR= HT/1000. 

C MLR = HM1/1000. 

BV( l) = 0. 

RV ( 2 ) = 1 . 

M B I M l = MB I - 1 

MB I P 1= MB I + 1 

MBOM 1= MBO-1 

M B 0 P 1 = MBO+1 

MMM1 = MM- l 

CP = AR/( GAM-1. 1*GAM 

EXPON= 1. /(GAM-1.) 

TWW= 2 . *OMEGA/ WTFL 

C PT I P= 2.*CP*TIP 

TGROG= 2.*GAM*AR/ l GAM + 1 • ) 

CALL SPLINT (MR, RMS P, NRSP , MV, MM, RM, SAL) 

CALI. SPLINT(MR,BESP,NRSP,MV,MM,BF,DBOM) 
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c 

C CALCULATE GEOMETRICAL CONSTANTS 
C 

CHORD! 2 ) = CHORD! 1) 

STGR (2 ) = STGR(l) 

ML E ( 1 ) = 0. 

MLE ( 2 ) = 0. 

THLF(l) = 0. 

THL E ( 2 ) = PITCH 
RMI ( 1) = RM ( MB I ) 

RMI ! 2 ) = R M ( M B I ) 

RMO ( I) = RM(MBO) 

RMO ( 2) = RM(MbO) 

C 

C INITIALIZE ARRAYS 
C 


CO 60 1=1,2500 
L ( I ) = 1. 

K I I ) = 0. 

60 RHO! I ) = RHOIP 

CO 70 IM* 1,100 
CO 70 SURF* 1 ,2 
RHOHB! IM, SURF) = RHOIP 
70 RHOVB! IM, SURF ) * RHOIP 
RETURN 

1000 FORMAT ( 1 H 1 ) 

1010 FORMAT ! 1615) 

1020 FORMAT (IX, 1617) 

1030 FORMAT (8F10.5) 

1040 FORMAT (1X.8G16.7) 

1100 FORMAT ( 80H 


1 

1110 FORMAT 
1 

1120 FORMAT 
1125 FORMAT 


) 

( 7X» 3HGAM, 14X, 2HAR, 13X , 3HT IP, 12X, 5NRH0IP, 12X,4HWTFL,11X,6H 
10X.5H0MEGA, 12X, 3H0RF) 

( 6X , 5HBETAI , 10X , 5HBE T At) , 1 IX , 6HCH0RDF, UX,5HSTGRF ) 

( 6X , 6HR EDFAC, 1 OX, 6H0ENT0L ) 


1130 FORMAT (41H MBI MQC MM NOB I NtiL NRSP) 

1140 FORMAT (39HL BLADE SURFACE l — UPPER SURFACE) 

1150 FORMAT ( 39HL BLADE SURFACE 2 — LOWER SURFACE) 

USO FORMAT ( 7 X , 2 HR [ , I 1 , 1 2 X , 2 HR 0 , II , 12X.4HBETI ,11,1 1X,4HBET0, I 1, 1 IX,5HS 
1 FLNG y 1 1 } 


1190 FORMAT ( 7X , 3HMSP , I 1, 2X , 5HARR AY ) 

1200 FORMAT !7X,4HTHSP, II , 2X, 5H ARRAY) 

1210 FORMAT (16HL MR ARRAY) 

1220 FORMAT (7X,llhRMSP ARRAY) 

1230 FORMAT ! 7X , 1 1HBESP ARRAY) 

1240 FORMAT (52HL BLOAT AANDK ERSOR STRTN SLCRD INTVL SURVL) 
1250 FORMAT (41H1 MM.NbBI , NRSP, OR SOMF SPLNO IS TOU LARGE) 

END 
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SUBROUTINE PR EC AL 


PRFCAL CALCULATES ALL REQUIRED FIXED CONSTANTS 


COMMON SRW»ITtR»IEND»LER(2),NER(2) 

COMMON / AUKRHO/ A ( 2500 , A ) , U { 2500 ) . K ( 2500 ) , RMO { 250U ) 

COMMON /INP/GAM, AR,TIP,RHOIP,WTFL,OMEGA,ORF,BLTAI .RtlAO.REUFAC, 

1 ^IJENTOL * Mb I » MBO » MM » N8B l , NPL » NRSP » MR ( 50 ) , RMSP (50) » R t SP ( 50) , 

2 KLD AT , AANDK t E^SGR v STR FN » SLC KD t I NTVL i SUR VL 

COMMON /CALCON/ACTWT, ACTGMG, *C TLAM.MB I M 1 ,M8 I P 1 , MBOMI , MB0P1 , M^M I, 

1 HMI t HT,DTLR,DMLR,PITCH.CP,EXPON,TVjW.CPTIP,TuROG,TBI,TdO,LAMBUA, 

2 TWL » l TM I N , I TM AX »N I P» I MS l 2 ) » BV ( 2 ) , MV ( 100 ) , I V 1 101 ) » I TV j 100 » 2 ) » 

3 TVaOO,2),UrOMvi l0 O f 2 I, HETAV1 100.21, MH1 100.2). OTOMH( 100,21. 

A 3ETAH( 100,2) »RMH( ICO. 2) ,8EH( 100,2) »RM( 100 ) ,BE( LOO) »OBDM( 10 » 

5 SALl 100) , AAA( 100) 

COMMON /02TDM2/ D2TOM2 ( 100 , 2 ) 

INTCGER 0 BLDAt!aANDK^ERSOR,STRFN,SLCRO,SURVL,AA TEMP, SURF, FIRST, 

PEAL K,KAK, LAMBDA, LMAX,MH,MLF,MR,MSL» MSP »MV, MV IMi 

EXTERNAL BLL.BL2 


CALCULATE LAMBDA, VI ANO VO 


C 

c 

c 

c 

c 


GET A I = BETAI/57.295779 

PET AO = BETAO/57. 295779 
TBI = S IN I BETA I ) /COS ( BET AI ) 

TBO = S I N ( 8c T AO) /COS ( BET AO ) 

10 RHOT = RHO I P , . 

RHOVI = WTFL/bE ( MB I ) /P ITCH/COS ( BETAI ) /RMIMBI ) 

20 VI = RHOVI/RHOT 

LAMBDA = RM { MB I ) • ( V I »SIN(BETAI ) +OMEGA*RM (MB I ) ) 

TTIP = l«~(Vl**2+2.»OMEG A* LAMBDA- ( OMEGA*RM ( MB I ) )**2) /CPTIP 
IF ( TTIP.LE.O. ) GO TO 30 
RHOMBI = RH0IP*TTIP**6XPCN 

I F ( ABS ( RHOMB I-RHOT )/ RFGI P. LT .. 000001 ) GO TO AO 
RHOT = RH0M3I 
GO TO 20 

30 WRI rE(6,102G) WTFL 
STOP 

AO VI = RHOVI /RHOMB I 

LAMBDA = RM ( Mb l) * ( VI * S IN (BET AI ) «-OMEGA*RM( MBI ) ) 

TWL = 2 . *OMtGA*L AMBD A 

RHOVO = WTFL/BUMBO ) / PI TCH/COS(BETAO)/RM(MBO) 

RH0 ! '*B2 = RHUIP 

TWLMR = TWL- ( OMEGA-*RM ( MBC ))**2 
LER ! I) = I 

CA f LL T JENSTY(pSovi,RHOMB2, VO, TWLMR, CPTIP, LXPON.RHOIP, GAM, AR, TIP) 

CALCULATE W-CRIT1CAL, AND MAXIMUM VALUES FOR RHO*W AT LEADING AND 
TRAILING EDGE 


AA ^ ( TWL-(UMEGA*RM( MBI ) )**2)/CPTIP 
TPP = T I P * ( I.-AA) 

PB = TGROG* TPP 
WCRI = SQRT(BB) 

TTIP = 1 . -BB/CPT I P-AA 

RHCrtMI = RHOIP»TTIP**EXPCN*WCRI 
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AA = (TWL-(OMEGA*RM(MeC CPTIP 

TPP = T I P » { 1 ,-AA ) ‘ 

PB = TGROG*TPP 
WCRG = SQRT(Bb) 

TTIP = l.-BB/CPTIP-AA 

PHOWMO = RHGIP*TTIP**EXPCN*WCRO 

vI"ES i o; U ;T F ri™ S L ^BSI FU a ' , ' :G '' 4N ° LW,,DA - “^ULAT,- REDUCED 

ACTWT = WTFL 
ACTCMG = OMEGA 
ACTLAM = LAMBDA 
WTFL = WTFL*REDF AC 
CMECA = OMEGA*REDFAC 


CALCULATE LAMBDA FOR REDUCED WEIGHT FLOW 


RHOT = RHOIP 

50 VT^= ^RHOvT/RHGT I/RMIMBI ) 

LAMBOA = RM(H8I).(VT*SIN(BETAI) + 0MEGA.RH(MBm 

RHO T ‘ RHOMB p »LT.. 000001. ) GO TO 60 

GO TO 50 

60 VT = RHOVI/KHGMBI 

LAMBDA = RM ( MB I ) * ( VT« SIN (BETAI )+OMEGA»RM IMB I ) ) 

TWL = 2.*0MEGA*LAMBDA 

S^GHt"^"* C0RRECTeC T0 8 °<™«RY »-« *N0 G-H USING REDUCED 
NERT = NER ( I ) 

TWL MR = TWL-(GMEGA*RM( 1 ) )*«2 
RH01 = RHOMB I 
TBI I = I.E20 

70 i IB ”wrfiIp!KH| l ™Ein H01/RH0M '* l * o " ,:G4 * ,R " ,MIJ,l " 2 - R, ' ,l, **^* ,i »o* 

IF( ABS ( TBI I-TB IT ) .LT . .00001) GU TO 80 
TBI I = TBIT 

LERaj=2 WTFL/PITCH#S0RT,1 ‘ +TBIU * 2, ' flEll,/RMll) 
censty CALL NO. 2 

GO L TO D 70 STV ,RHC1VI »RHCl, AA,TWLMR f CPTIP, EXPON.RHO IP, GAM, AR, TIP) 

80 TBI = TBIT 

Bwnv^ = - W I FL/BCU ' Bn ,/PITCH /COS(BETAO)/RM(MBO) 

KnUi o 2 — RHOIP 

TWLMR = TWL-(OMEGA*RM(MBC ) ) «*2 
LER ( 1 ) = 3 

CENSTY CALL NO. 3 

efiL^mml[?:5y?;S5-;“ ,T " LKR,c " l,, ' EXI ’ ON ' iiHoi ' , * i: *">*«-Tipi 

TWL MR = TWL-(0MtGA*RM(MM))**2 
RHOMM = RHOM02 
T BOR = 1.E20 

VO TBGT - (TBO/BEIMBO ) * RH0MM/RH0MB2+0MEGA* (RM( MBO J **2-RM ( MM ) **2 ) * 
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1 PHOMM/WTFL*PITCH)*BE(J^M) 

TBo”° S ™? M ~ r80T, ' Lr " OCOO ‘’ G0 Tn 100 

CENSTY CALL NU. 4 

cS'VrT"' IRH0V0 - RHD ""' 4 ».™L"R.CPTIP.txP0N,KH0IP,GAP,AR.TIP| 

100 TBO = TBOT 

BTAUUT = ATAN ( TBfJ) *5 7 .295779 
I F ( N E R ( D.tQ.NERT) GO TO 110 
hRIfE (6,1025) 

STOP 


CALCULATE TV, [TV, IV. DTDM V, AND BETAV ARRAYS 


110 ITMIN = 0 

ITMAX = NOBI-i 

C TV, I TV, AND DTDMV ON BLADE 
CO 120 IM=M6I,M80 
LER ( 2) = 1 

C PLCD CALL NO. I 

CALL BL 1 ( M V ( IM ) , TV ( I M , l ) ,DTOMV (IM , 1 ) . INF » 

TVHM.l), INT((TV(I M ,1) + DTLR)/HT) 

IF ( TV f IM, 1 ) .GT.-OTLR ) IlVdM. 1) = ITVI IM 11*1 

ITMIN= MINOIITMIN.ITVlIH.l ITVdM.im 

L ER ( 2 ) =2 

C BLCD CALL NO. 2 

^™!\m L o! MV(IM),TV<IM ’ 2, ’ DT0MV (I m *2>.INF) 

ITV(IM,2)« INTI (TV(IM,2)-DTLR)/HT) 

IE (TV(IMf2)*LT.0TLR) ITV(IM*2)=!TVfrM ? ) i 
r 120 I TM AX= MAX0dT.MAX,ITV(IM,2n ' ’ 2,_l 

- I TV AND IV UPSTREAM OF eLACE 
FIRST = 0 
LAST = N8B I- 1 
CO 130 I M= l , Mb I M 1 
ITVl IM, 1)= FIRST 
130 I TV ( IM, 2) s LAST 
: ITV DOWNSTREAM OF BLADE 

1^0 LAST= I TV ( MdO , 2 ) 

FIRST= IAST+1-NBBI 
CO 150 IM=MBOP 1 , MM 

I TV ( IM , 1 ) = first 
150 I TV ( I M , 2 ) = LAST 

MINOdTMIN.ITVCMM, 1) ) 

CALCULATE IV ARRAY 
IVd ) = l 
CO 160 IM= 1 , MM 

CO 200 SURF= 1 , 2 
CO 200 IM = Mci I , MBO 

l cuK n:::^^ i i,: 0 ^;7^;?ss^f“5;^rr SAL, ' ,M ’* DTo " v<iM ' suRF>> > 

M"' , :i;&,T N,OTn " vun>SURF, * RB ' , " i) * 57 -” G ^ 

p|;"™“ R f I | ;j" , - eT4 ' N ''' n ' , < H <>“™.-CRO,8TAaor 

, K o 5 If 5 1 tl°50) ITMIN, I TMAX.ACTLAM, LAMBDA, NIP 
WR I TE ( 6 , 1 060 ) ( SURF, BV (SURF) ,SURF=I,2J 
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IF(BLDAT.LE.O) GO TO 2 30 

WRITE ( ft! 1060 ) (MVllM),TV!IM,l> » DTDMV ( IM, I ) > CURVI IM,i),TV(IM,2>, 

1 OTDMV l I M , 2 ) t CURVl I M • 2 ) , IM=MR I . MBO ) , uu » 

WRITE 16, 109 0) ( IM,MV!IM),RM( IM) ,SAL( IM) ,BF! IM ) , DOOM ( I M) , IM-l.MM) 

230 CONTINUE 


C 

c 

c 


CALCULATE MH AND DTDMH ARRAYS 


ITO = ITVU.l) 

MRTS = 1 
IMS ( 1 ) = 1 
MH ! 1 » 1 ) - 0. 

CTDMH! 1,1) = 1.E10 
LERI") = 3 

BLCO AND ROOT (VIA MHCRIZ) CALL NO. 3 ull ,, 

CALL MHORIZIMV, I TV ( 1 , 1) , ELI, MB I, MBO, I TO , HT ,DTLR , 0 , I MS ( 1 ) ,MH( l, 1) , 

1 DTDMH ( 1,1), MRT S ) 

IF I ITVIMBO, 1 )-ITV(MR0,2 1+NBBI .NE.2 ) GO TO 240 

I MSL = IMS! 1 ) + l 

MH ( I MSL , 1 ) = MV I MBO ) 

CTDMHI IMSL.i) = -1.E10 
I MS 1 1 ) = IMSL 
240 IMS12) = 0 
MRTS = 1 
LERI2) =4 

BLCD AND ROOT (VIA MHORIZ) CALL NO. 4 MMM ,, 

CALL MHORIZ (MV, IT V(l, 2) , BL2 , MB I , MBO, I TO , HT , DTLR , 1 , I MS ( 2 ) , MH ( l , 2 , 

1 DTDMH( 1,2) »MRTS) 

I = MAXO ( I MS ( 1 ) • I M S ( 2 ) ) 

IF ( I .LE. 100 ) GO TO 29C 

WRI TE( 6, 1 100) I 

STOP 

290 WRI TE°t 1 110) ° UM° I V( IM) , ( ITV( I M , SURF ) , SURF= 1 , 2 > , IM=1,MM) 


CALCULATE RMH , BEH* AND BET AH ARRAYS 


300 IF (BLDAT.GT.O) WR I TE ( 6 , 1 120 ) 

CALL SPLINT! MR, RMSP, NR SP,MH( 1, SURF ) , IMS ( SURF ) »RMH< l , SURF ) »**A) 
CALL SPLINT! MR, BESP,NRSP,MH( 1 , SURF ), IMS 1 SURF >, BEH ( l , SURF ) ,AAA) 
IMSS = IMS(SURF) 

IF ( IMSS.LT. 1 ) GO TO 320 

310 eETAH! I Hs! SURF ) '= S AT AM DTDMH UHS, SURF) *RMH! I HS, SURF) ) *57.295779 

5f JSlDAtIgT.O) WRITE 16,1130) SURF ,( MH! IM, SURF ), RMH! IM. SURF ) , 
l BEH ( I M, SURF ) ,BETAH l IM, SURF) .DTDMH! IM, SURF ), IM=1 , IMSS) 

320 CONTINUE 

IF (BLDAT.Lt.O) GO TO 340 
WRITE (6,1140) 

IT = ITMIN 

330 IF (IT.GT.ITMAX) GO TO 340 
TH = FLOAT! IT)*HT 
WRITE (6,1010) I T , TH 
IT = IT+l 
GO TO 330 

340 IFINIP.LE.2500) GO TO 350 
WRITE(6,1150) 

STOP 
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350 in R I T E 16,1000) 

RETURN 

1000 FORMAT (1H1) 

1010 FORMAT ( 4X , 1 4 , G 1 6 . 5) 

1020 FORMAT l 6 OHL I N PUT WEIGHT FLOW ( WTFL ) IS TOO LARGE AT BLADE LFADING 
1 EDGE ) 

1025 FORMAT ( 40HL REDUCED WEIGHT FLOW IS STILL TOO LARGE) 

1030 FORMAT { IHI/24X, 10HF R EES TR E AM , 8X , 13HMAXIMUM VALUE, 

17X,HHCRITICAL,30X, 14 H BETA CORRECTED/ 25 X, 8HVEL0C ITY»10X»9HF0R RHO*W 
2, 10X.3HVELOC ITY, 31X, 1 1HT0 BOUNDAR Y/ 1 X , 1 7HLE AD I NG EDGE B-G,3Gl6.5, 
312X, 12 H BOUNDARY A-H, G 1 8 . 5/ IX , 1 7HTRA I L ING EDGE C-F , 3G 1 B . 5 , 12X , 

4 1 2 H rt OU N D A R Y D-E , G1 8. 5 / 86 X, 30H ( BASED ON REDUCED WEIGHT FLOW)) 

1040 FORMAT ( 33HL CALCULATED PROGRAM CONSTANT S//5X,5HPITCH, 13X, 

1 ?HHT,13X,3HHMl/lX,5G16.7) 

1050 FORMAT ( / 5X , 5H I TM l N, l OX , 5H I T MAX/ 4X , I 5 , 10X, I 5 // 5X , 6HL AMBD A , 12X, 

1 29HL AMBOA AT REDUCED WEIGHT FLOW/ IX , G 16 . 7 , 1 2X ,G 1 6. 7/ 

2 38HL NUMBER OF INTERIOR MESH POINTS = ,13) 

1060 FORMAT ( 28HL SURFACE BOUNDARY VALUES// 5X, 7HSURF ACE , 7X , 2HBV 
1/(5X,I4,4X»F10.5) ) 

1070 FORMAT ( 1 H 1 , 6X , 62 HBL AGE DATA AT INTERSECTIONS UF VERTICAL MESH LIN 
IES WITH BLADFS) 

1080 FORMAT ( IHL , 22X, I5HBL ADE SURFACE 1 , 30X , l 5HBLADE SURFACE 2/7X, 

1 1HM, 14X.2HTV, l IX, 5 HOT DM V, l IX , 4HCURV, 1 2X , 2HTV , 1 1 X , 5HDTDMV , 1 1 X, 

2 4HCURV/I7G15.5) ) 

1090 FORMAT ( 1 H I , 1 3X , 44HSTR EAM SHEET COORDINATES AND THICKNESS TABLE / 

1 2X,2HIM,7X, IHM, 14X.1HR, I3X, 3HSAL , 13X, IHB, 12X, 5HDB/DM/ ( LX, 13, 

2 5G15.5)) 

1100 FORMAT I 34HL0NE OF THE MH ARRAYS IS TOO LARGE /7H IT HAS, 15, 8H POI 
1NTS ) 

1110 FORMAT (4H1 IM,9X,8HIV ARR AY , 25X , 9H I TV ARR A Y/38X , 5HBL ADE/ 37X , 7HSUR 
IFACE.3X, 1H1 , 5X , IH2/3 9 X , 3HN0. / ( IX , 13, 5X , I 10 , 25X , 2 ( I 4, 2X ) ) ) 

1120 FORMAT (67H1M COORDINATES OF INTERSECTIONS OF HORIZONTAL MESH LINE 
IS WITH BLADE) 

1130 FORMAT (25HLMH ARRAY - BLADE SURFACE, I2//15X,2HMH,19X,3HRMH, 19X, 
l 3HBEH » 1 8X » 5HBET AH » 1 7X , 5HDT0MH/ ( 5G22 .4 ) ) 

1140 FORMAT (43 H IT H ETA COORDINATES OF HORIZONTAL MESH LINES// 6X.2HIT, 

1 5X , 5HTHET A ) 

1150 FORMAT 1 48HLTHE NUMBER OF INTFRIOR MESH POINTS EXCEEDS 2500) 

END 


SUBROUTINE COEF 

COE F CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K, 
AT ALL UNKNOWN MESH POINTS FOR THE ENTIRE REGION 

COMMON SRW , ITER, I END, LFR (2 ) , NER( 2) 

COMMON / AUKAHO/ A ( 2500 , 4 ) , U ( 2500 ) , K ( 2500 ) , RHO ( 2500 ) 

COMMON /INP/GAM, AR,T I P , RHO I P , WTF L , OMEGA , ORF , BE TA I , BE T AO, REDF AC , 

1 DENTOL ,MB I , MBO , MM , NRB I , NB L , NRSP , MR ( 50 ) , RMSP l 50 ) , BESP ( 50 ) , 

2 BLDAT , AANDK , ERSOR , STRFN, SLCRD, INTVL.SUAVL 

COMMON /CALCON/ACTWT, ACT CMG, AC TL AM , MB I Ml ,MH I P 1 , MBOM I , MBO PI , MMM 1 , 

1 HM1 , HT,DTLR,OMLR, PITCH, CP, EXPO N,TWW,CPTIP»T GROG , TB I ,TBO, LAMBDA, 

2 TWL, ITMIN, ITMAX.NI P, IMS (2 ) , 8V ( 2 ) , MV ( 1 00 ) , I V ( 10 1 ) , I T V ( 100 , 2 ) , 

3 TV( 100 , 2 ) , LT DMV (100, 2), BET AV (100, 2), MH 1100, 2) , DTOMH( 100,2 ) . 

4 BET AH ( 100,2) , RMH( 100, 2) ,BFH( 100 , 2 ) , RM ( 100 » , BE ( 1 00 ) , DBDM ( 100 ) , 

5 SAL( 100) ,AAA( 100) 


COMMON /HRBAAK/HI4), R (4) ,B (4) ,KAK(4) , KA(4) »RZ»bZ»IH(4) 

INTEGER BLOAT, AANDK, FRSOR, STRFN, SLCRD, SURVL , AA TEMP, SURF, FIRST, 
1 UPPER, SI, ST, SRW 

REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1 
C INITIALIZE ARRAYS 
ITER = I TER-*- 1 
IH( l } = 1 

IH{2) = 0 

C INCOMPRESSIBLE CASE 

IF (GAM.NE. 1.5.0R. AR.NE. 1000. .OR.TIP.NE. 1.E6) GU TO 20 
I END = 1 
GO TO 40 

C ADJUSTMENT OF PRINTING CONTROL VARIABLES 
20 IF( ITER.NE. I. AND. ITER.NE .2) GO TO 30 
AANDK = AANDK- l 
ERSOR = ERSOR-1 
STRFN = STRFN- 1 
SLCRD = SLCRD- 1 
INTVL = INTVL-1 
SURVL = SURVL-1 
30 IF ( I END.NE .0 ) GO TO 40 
AANDK = AANOK+2 
ERSUR = ERSOR+2 
STRFN = STREN+2 
SLCRD = SLCRD+2 
INTVL = INTVL+2 
SURVL = SURVL+2 
C 

C FIRST VERTICAL MESH LINE 
C 

40 CO 50 I P= 1 , NBC I 
A ( IP, I) = 0. 

A ( I P , 2 ) = 0. 

A ( IP, 3 ) = 0. 

A ( I P ,4 ) = 1. 

50 K(IP) = HMI*TBI/PITCH/( CRM(I)+RM(2) J/2. ) 

C UPSTREAM OE BLADE, EXCEPT FOR FIRST VERTICAL MESH LINE 

IF (2.GT.MBIM1 ) GO TO 70 
CO 60 IM=2 ,MB I Ml 
60 CALL COEFPI I M ) 

C 

C BETWEEN BLADES 
C 

70 CO 30 I M= MB I , MBO 
80 CALL COEFBB ( IM) 

C 

C DOWNSTREAM OF BLADES EXCEPT FOR FINAL MESH LINE 
C 

150 IFIMB0P1.GT.MMM1 ) GO TO 170 
CO 160 I M=MbOP 1 , MMM1 
160 CALL COEFPI IM) 

C 

C FINAL VERTICAL MESH LINE 
C 

170 I VMM = I V ( MM ) 

CO 180 I P= I VMM , N I P 
A ( I P , 1 ) = 0. 

A( IP, 2) = 0. 
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A ( I P , 3 ) = 1. 

A( IP, 4) = 0. 

180 K ( I P ) = -HMl*TBO/PITCH/RM(MM) 

TAKE CARE OF POINTS ADJACENT TO B, AND CASES WHEN POINTS J.C,E, OR F 
ARE GRID PniNTS 

POINT B 

IP = IV(MBIMl) 

A ( IP, 4) = 0. 

POINT C 

IF( ITV (MBO. 1 )-ITV(MBO t 2) +NBBI.NE.2) RETURN 
IT = I TV ( MBQ » 1 ) — 1 
IP = I PF ( MBOPL , IT) 

A ( IP, 3) = 0. 

RETURN 

END 


SUBROUTINE COEFBB ( I M ) 

COEFBB CALCULATES FINITE DIFFERENCE COEFFICIENTS. A, AND CONSTANTS, K 
ALONG ALL VERTICAL MESH LINES WHICH INTERSECT BLADES 

COMMON / AOKRHO/ A ( 25 00 , A ) , U ( 2500 ) , K ( 2500 ) , RHO I 2500 ) 

COMMON /INP/GAM, AR,T I P.RHOIP ,WTFL» OMEGA ,ORF , BE T A I » BETAO , REDF AC . 

1 DENT OL, Mo I, MBO , MM , NBB I , NBL , NRSP , MR ( 50 J.RMSPI50) » BESP I 50 ) , 

2 BLDAT . AANDK , ERSOR « STRFN. SLCRD. I NTVL. SURVL 
COMMON /C ALCON/ACTWT , ACT CMG , AC TL AM , MO IM 1 .MB I P 1 , MBOM 1 , MB0P1 . MMM1 , 

1 HM1 .HT.DTLK.DMLR, P I TC H, CP , E XPON , TWW , CPT I P .TGROG , TB I ,TBO. LAMBDA, 

2 TWL » ITMIN, 1TMAX.NIP, I MS (2 ) »BV( 2 ) »MVl 100) , IVI 101 ) , IT VI 100,2) , 

3 TV! 100 , 2 ) , DTOMV ( 100, 2 ) , BETAV l 100, 2 ) , MH( 100,2) , DTD MH (100,2) , 
if BETAHl 100,2) , RMH( 100,2) ,BEH(100,2)»RM(100),BE(100) , OB DM (100), 

5 SAL( 100) , AAAI 100) 

COMMON /HRRAAK/H(4),R(4) ,B ( A ) , KAK ( A) ,KA(4J,RZ,BZ,IH(4) 

INThGER BLDAT, AANDK, ERSOR, STRFN , SLCRD , SURVL , AATEMP , SURF , F I R S T , 

1 UPPER, SI, ST, SRW 

REAL K,KAK,LAMBDA,LMAX,MH,MLE.,MR,MSL,MSP,MV,MVIM1 

I F ( I TV l IM, 1) .GT.ITVt IM,2 )) RETURN 
(TVU = I T V ( I M , I ) 

IT = I TVU - 1 
ITVL = I TV ( I M , 2 ) 

IPU = I PF ( [M, ITVU) 

I PL = I PU+ I TVL-I TVU 
CO 90 I P= I PU , I PL 
IT = I T + 1 

CALL HRO ( IM, IT , IP) 

CO 10 1=1,4 
K AK ( I ) = 0. 

10 K A ( I ) = 0 

C FIX HRB VALUES FOR CASES WHERE MESH LINES INTERSECT BLADES 
60 IF! IT.EQ. ITV( IM, 1) ) CALL BDR Y 1 2 ( 1 , I M, I T ) 

I F ( IT.EQ. I TV ( I M , 2 ) ) CALL BDR Y 1 2 ( 2 , I M , I T ) 

ITVM1 = I TV ( IM-1 , 1 ) 

ITVPl = I TV ( I M+ 1 » 1 ) 

I F ( IT.LT. ITVM1) CALL BDR Y34 ( 3 , I M , 1 ) 
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IF! IT.LT. ITVP1 ) CALL BDR Y34 ( 4, IM, 1 ) 

IF( IT.GT. ITVIIM-1,2) ) CALL BDRY34 I 3, IM, 2 ) 

IFIIT.GT. ITVIIM+1,2) ) CALL BDRY34 l 4, IM, 2 ) 

70 I F { IM.EQ. MBO. AND. LOWER . EC. 2 ) GU TO 80 
C COMPUTE A AND K COEFFICIENTS 
80 CALL A AM IM, IP) 

CO 90 1=1,4 

K I I P ) = K(IP)+KAK(I)*A(IP,I) 

90 IF(KAm.EQ.l) A ( I P , I ) * 0 . 

RETURN 

C 

C COEFP CALCULATES FINITE C I F FERENCE CUEFFI C I ENTS , A, AND CONSTANTS, 
C ALONG ALL VERTICAL MESH LINES WHICH DO NOT INTERSECT BLADES 

ENTRY COE FP l IM) 

ITVU = ITVI IM, l) 

IT = ITVU- l 
IT VL = I T V ( I M , 2 ) 

I PL = I V { IM+1J-1 
IPU = I V ( IM) 

CG 100 I P = IPU , I PL 
IT = IT+1 

CALL HRB ( IM, IT, IP) 

IF UT.EQ.ITVU) R I 1 ) = RHO(IPL) 

IF UT.EQ.ITVL) R C 2 ) = RHOIIPU) 

100 CALL A AK ( IM, IP ) 

KC I PL) = K ( I PL ) + A { I P L , 2 ) 

K I IPU) = K ( I PU ) -A ( IPU , 1 ) 

RETURN 

END 


K, 


SUBROUTINE HRB I IM, IT, IP) 

HRB CALCULATES MESH SPACING, H, DENSITIES, RZ AND R, AT GIVEN AND 
ADJACENT POINTS, AND STREAM SHEET THICKNESSES, BZ AND B, AT GIVEN 
AND ADJACENT POINTS 

COMMON / AUKRHG/ A ( 25 00 , 4 ) , U ( 2500 ) , K ( 2500 ) , RHO ( 2500 ) 

COMMON /CALCON/ACTWT , ACTOMG, AC TEAM, MB I Ml ,MBIP1 ,MBOMI ,MBOPi , MMM I , 

1 HMi ,HT,DTLR» DMLR, P I TC H, CP , EXPON ,TWK,CPTIP,T GROG ,TBI,T BO, LAMBDA, 

2 TWL, ITMIN,ITMAX»NIP,IMS(2),8V(2)»MV( I00),IV(10l),ITV( 100,2) , 

3 TVl 100,2) ,DTOMV( 100,2 l.BETAVI 100,2) ,MH( 100, 2 ) . DTDMH { I 00 , 2 ) , 

4 BET AH ( 100,2 ) ,RMHl 100, 2) , BEH { 100,2) »RM( 100) »GE ( 100) ,DBDM( 100 ) , 

5 SAL ( 100) , AAAI 100) 

COMMON /HRBAAK/HI 4 ) , R ( 4 ) ,B l 4 ) , K AK ( 4 ) ,KA ( 4) , RZ , BZ , I H { 4 ) 

INTEGER BLDAT, AANDK, ERSOR, STRFN, SLCRD, SURVL , AATEMP , SURF , FIRS T, 

1 UPPER, SI, ST, SRW 

REAL K,KAK,L AMBDA , LM AX , MH, ML E , MR , M SL , MSP , MV, MV IM 1 
HI 1 )= HT* RM I IM ) 

H l 2 ) = HT*RM ( IM) 

H I 3 ) = MV I IM) - MV I IM- l ) 

H ( 4 ) = MVI IM+i )-MV( IM ) 

RZ = RHOI IP) 

IP3 = IPFI I M - 1 , IT) 

IP4 = IPFI IM + 1, IT) 
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R ( L ) = RHOIIP-1) 
R(2» = RHO ( I P •*- 1 ) 
R { 3 ) = RHO(IPi) 

R ( 4 ) = RHO I I P4 ) 
BZ= BE ( I M ) 

B ( 3 ) = BE ( I M- 1 ) 

E (4 )= BE ( IM+1 ) 

RETURN 

END 


SUBROUTINE AAKIIM.IP) 

AAK CALCULATES FINITE DIFFERENCE COEFFICIENTS. A, AND CONSTANT, K, 

AT A SINGLE MESH POINT 

COMMON / AUKRHQ/ A ( 250C , 4 ) , U I 2500 I , K ( 2500 ) , RHO 1 2500 ) 

COMMON /CALCON/ACTWT, ACTOMG, ACTLAM.MBIMl .MB I PI .MB0M1 .MBOPI » MMM1, 

1 HM1 .HT.DTLR.DMLR. P ITCH, CP . EXPON , TWW, CPT I P , TGROG » TB I . T BO, LAMBDA, 

2 TWL, ITMIN.ITMAX.NIP, IMS12 ) ♦ BV 1 2 ) , MV I LOO ) , I V 1 101 ) , I TV! 100 , 2 ) , 

3 TV( 100,2) ,UTDMV(100,2),BETAVl 100,2) ,MH( 100 , 2 ) , DTDMH ( 1 00 , 2 ) , 

4 BET AH ( 100,2 ) , RMHI 1 00 , 2 ) , BEH ( 100 , 2 ) , RM I 100 J , t»E I 1 00 ) , DBDM ( 100 ) , 

5 SAH 100) , AAAI 100) 

COMMON /HRBAAK/H(4),R (4) ,BI4) ,KAK14) ,KA(4) , RZ , BZ » IHI4) 

INTEGER BLOAT, AANDK, ERSOR, STRFN, SLCRD, SURVL , AATEMP , SURF , F I RST , 

1 UPPER, SI, ST, SRW 

REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP.MV,MV1M1 
A 12 - 2. /HI D/H12) 

A34= 2./HI 3)/HI4) 

AZ= A12+A34 

B 1 2= (R12)-RI1))/RZ/IHI1)+H(2) ) 

034= (BI4)*R(4)-BI3)«R(3))/6Z/RZ/IH(3)+H(4) J-SALI IM)/RMI IM) 

A I I P , 1 ) = 12. /HI 1 )+B12l/AZ/(H( 1)+HI2) ) 

A I I P , 2 ) = A12/AZ-AI I P , 1 ) 

AIIP.3) = I 2. /HI 3)+B34 ) / AZ/ IHI 3 ) + H I 4) ) 

A ( I P,4 ) = A34/AZ-AI IP, 3) 

K ( I P ) = -TWK*BZ*RZ*SALIIM)/AZ 

RETURN 

FND 


SUBROUTINE 6DRY12 I I, I M , I T) 

C 

C BDRY12 CORRECTS VALUES COMPUTED BY HRB WHEN A VERTICAL MESH LINE 
C INTERSECTS A BLADE 
C 

COMMON /CALCON/ACTWT, ACTOMG, ACTLAM.MBIML, MPIP1.MB0M1 , MBOPI, MMM1, 

1 HM1, HT.DTLR.DMLR, PITCH, CP, EXPON, TWW , CPT I P, TbROG , TB I , T BO, LAMBDA, 

2 TWL, ITMIN.ITMAX.NIP, I MS I 2 ) , BV I 2 ) . MV I 100 ) , I V 1 1 0 1 ) , I I V ( 1 00 , 2 ) , 

3 TV! 100,2) , DTOMV 1 100,2 ) , BETA VI 100, 2 ) ,MH( 100, 2), DTDMH 1100, 2) , 

4 BET AH I 100, 2), RMHI 100,2), BEHI 100,2), RMI1 00 ),BEl 100) ,DBDMI 100), 

5 SAL! 100) , AAAI 100) 

COMMON /RHOS/RHOHBI 100, 2 ),RH0V3I 100,2) 
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COMMON /HRBAAK/H(4),R(4) ,B { 4 ) , KAK ( 4 ) f KA ( 4 > ,RZ , 8Z , IH ( 4 ) 

INTEGER BLDAT, AANDK, ERSOR, STRFN »SLCRD« SURVL , A A TEMP , SURF ,F I RST, 
I UPPER, SI, ST, SRW 

REAL K, KAK, LAMBDA, LMAX,MH,MLF,MR,MSL, MSP, MV , MV IM I 
F(I) = ABS ( FLOAT ( I T ) *HT-TV (IM, I ) ) *RM( IM ) 

R( 11= RHOVB ( IM, I ) 

KAK ( I ) = 8 V { I ) 

K A ( I ) = 1 

RETURN 

END 


SUBROUTINE B0KY34 ( I , I M , SURF ) 

C 

C BORY34 CORRECTS VALUES COMPUTED BY HRB WHEN A HORIZONTAL MESH LINE 
C INTERSECTS A BLADE 
C 

COMMON /C ALCON/ACTWT, ACTCMG, AC TL AM, MB I MI, MB I Pi ,MBOMI , MBOP1 , MMM I , 

1 HMI , HT , DTLR , DMLR, P I TCH, CP , EXPON , TWW , CPT I P , TGROG , TB I ,TBO, LAMBDA, 

2 r WL, ITMIN»ITMAX,NIP,IM$(2)»6V(2)»MV( 100 ), IV( 101 ) , ITV( 100 , 2 ) , 

3 TV ( 100,2) ,DTDMV ( 1 00 , 2 ) , BETAV { 1 00 , 2 ) , MH( 100, 2 ) , DTDMH { 100,2), 

4 BET AH ( 100.2) ,RMH< 100, 2) ,BEH( 100,2) ,RM( 100) ,bE ( 1 00 ) , DBDM { 100 ) , 

5 SAL( 100) ,AAA( 100) 

COMMON /RHOS/RHOHBI 100,2) » RHOVB I 100,2) 

COMMON /HRBAAK/H(4),R(4) ,B(4),KAK(4),KA(4) f RZ,BZ,IH(4) 

INTEGER BLDAT, AANDK, ERSOR, STRFN, SLCRD, SURVL , A A TEMP , SURF , F I RST , 

I UPPER, SI, ST, SRW 

PEAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL .MSP ,MV, MV I M I 
I H { SURF )= IH( SURF ) +1 
I HS= IH ( SURF ) 

F ( I ) =ABS ( MV ( I M ) -MH ( IHS.SURF) ) 

R( I ) =RHOHB ( IHS.SURF) 

B ( I ) =BEH [ IHS.SURF ) 

KAK(I) = BV(SURF) 

KA ( I ) = 1 

RETURN 

END 


SUBROUTINE SOR 

SOR SOLVES THE SET OF SIMULTANEOUS EQUATIONS FOR THE STREAM FUNCTION 
USING THE METHOD OF SUCCESSIVE OVER-RELAXATION 

COMMON / AUK RHO/ A ( 25 0 0 , 4 ) , U I 2500 ) , K ( 2500 ) , RHO ( 2500 ) 

COMMON / I NP/GAM, AR , T I P ,RhO IP , WTFL , OMEGA ,ORF , BETA I , BETAO, REDF AC , 

1 DENTOL , MB I , M30,MM ,NBBI, NBL, NRSP»MR( 50 ) »RMSP ( 50) ,BESP( 50) , 

2 BLDAT, AANDK, ERSOR, STRFN, SLCRD, IN TVL, SURVL 

COMMON /C ALCON/ACTWT ,ACTCMG,ACTLAM,MBIMI,MBIP1 ,MROMI,MBOPI ,MMM1 , 

1 HMI, HT, OTLR, D ML R, PITCH, CP, EXP ON, TWW.CPTIP, TGROG, TBI,T BO, LAMBDA, 

2 TWL , ITMIN, l T M A X » N I P , I MS { 2 J , BV ( 2 ) , MV ( 100 ) , I V ( 10 1 ) , I TVI 100, 2 ) , 

3 TV { 100,2) .DTDMV l 100,2 ), BETAV ( 100, 2),MH( 100, 2 ) , DTDMH ( 100 , 2 ) , 

4 BET AH( 100 ,2 ) , RMH( ICO, 2) ,BEH( 100,2 ) ,RM( 100) ,BE { 100) ,OBDM ( 100 ) , 
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5 SAL ( 100) ,AAA( LOO) 

INTEGER BLDAT, AANDK, ERSOR, STRFN , SLCRD , SURVL . AATEMP , SURF , F I RST , 

1 UPPER, SI, ST, SRW 

PEAL K,KAK, LAMBDA, LMAX,MH,MLE, MR, MSL, MSP »MV , MV I Ml 
A AT EMP = AANDK 
IF (0RF.GE.2.) ORF=0. 

IFIORF.GT. 1. ) GO TO 50 
CRF = 1. 

CRFOPT = 2. 

AO CRF I EM=ORFOPT 
LMAX = 0. 

50 IF ( AATEMP.GT.O ) WR IT E ( 6 , 10 1 0 ) 

ERROR = 0. 

SOLVE MATRIX EQUATION BY SOR, OR CALCULATE OPTIMUM OVERRELAXATION 
FACTOR 

IP = 0 

CO 120 l M = 1 , MM 
IPU = I V ( I M ) 

I PL = I VI IM+l )-l 
IT * I T V ( IM, l) 

IF ( AATEMP.GT.O ) WRITE (6,1020) IM,IT 
CO 120 I P= I PU , IPL 
I PI = I P-1 
IP2 = IP+1 

C CORRECT IP1 AND I P2 ALONG PERIODIC BOUNDARIES 
I F C IM.GE.MBI .AND. IM.LE.MBO) GO TO 60 
IF! IT.EQ. ITVIIM, 1 ) ) IPl = IPl+NBBI 
IF! IT.EQ. I TV ( 1M,2> ) IP2 = IP2-NB8I 
60 IT3 = IT 
ITA = IT 

100 IP3 = I PF t IM-1, I T 3 ) 

IPA = I PF C IM+l, ITA) 

IFIORF.GT. 1. ) GO TO 110 
C CALCULATE NEW ESTIMATE FOR LMAX 

UNEW = A(IP,1)*U(IP1)+A(IP,2)*U(IP2)+A(IP,3)*U(IP3)+AUP,A)*U(IPA) 

IF (UNEW.LT.l.E-25) U ( T P ) = 0. 

IF (U( IP) .FU.O. ) GO TO 115 
RATIO = UNEw/Ul IP) 

LMAX= AMAX1 ( RAT 10, LMAX ) 

UUP) = UNEW 
GO 10 115 

C CALCULATE NEW ESTIMATE FOR STREAM FUNCTION BY SOR 

110 CHANGE = ORF*(K( IP )-U ( IP > + A( IP, 1 )*U( IPl )+A{ IP,2)*UI I P2 )+A( IP, 3)« 

1 U ( I P 3 ) + A ( I P , A ) *U ( IPA) ) 

ERRUR= AMAX1(ERR0R,ABS(CHANGE) ) 

L( IP) = U( I P ) +CH ANGE 
115 IF(AATEMP.LE.O) GO TO 120 

WRITE (6, 1030) I T, IP, IPl ,IP2,IP3, IPA, (A( IP, I ) , 1 = 1, A) ,K( IP) 

120 IT = IT+1 
AATEMP = 0 

IFIORF.GT. 1. ) CO TO 130 

CRFGPT= ? ./{ 1 . +SQRT ( A BS ( 1 . -L MAX ) ) ) 

WRITE(6,1000) ORFOPT 

IF (ORFTEM-ORFOPT.GT. .00001. OR. ORFOPT.GT. 1.999) GO TO AO 
WRITE (6,1070) 

CRF = ORFOPT 
GO TO 50 
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130 IF (ERSOR. GT.O ) WR I TE { 6 , 1 040 ) ERROR 
IF( ERROR. GT. .000001) GO TO 50 
IF ( STRFN. LE.O ) RETURN 
C 

C PRIM STREAM FUNCTION VALUES FOR THIS ITERATION 
C 

WRITE (6,1050) 

CO 140 IM=1,MM 
IPU = IV( IM) 

I PL = I V ( I M+ 1 ) — 1 
ITVU = IT V ( IM, l) 

WRITE (6,1020) I M , ITVU 
140 WRITE ( 6 , I0t>0 ) ( U ( IP ) , I P = I PU , I PL ) 

RETURN 

1000 FORMAT ( 24H ESTIMATED OPTIMUM ORE =,F9.6) 

1010 FORMAT (82H1 IT IP I P 1 IP2 IP3 IP4 All) A(2) 

1 A{ 3 ) A ( 4 ) K) 

1020 FORMAT ( 5HK IM = , I 4 , 6X , 6H I T1 = ,14) 

1030 FORMAT ( IX , 14 ,5I6,5F10 .5) 

1040 FORMAT ( 8H ERROR =,F11.8) 

1050 FORMAT ( 1H1, 10X.46HSTREAM FUNCTION VALUES FOR REDUCFD WEIGHT FLOW) 
1060 FORMAT (2X, 10F13.8) 

1070 FORMAT (1H1) 

END 


SUBROUTINE SLAX 
C 

C SLAX CALLS SUBROUTINES TO CALCULATE RHG*W-SUB-M THROUGHOUT THE REGION 
C AND ON THE BLADE SURFACES, AND TO CALCULATE AND PLOT THE 
C STREAMLINE LOCATIONS 
C 

COMMON / AUKRHO/ A ( 2500 , 4 ) , U I 2500 ) , K ( 2500 ) , RHO( 2500 ) 

COMMON /INP/GAM,AR,T IP , RHO I P , WTFL , OMEGA , ORF , BE TA I , BF T AO , REDF AC , 

1 DENT OL , MB I » MBO»MM » NBB I » NBL , NRSP, MR { 50 > » RMSP ( 50 ) , BESP { 50 ) » 

7 DLDAT , A ANDK , ERSOR , S TR EN , SLCRD , I NTVL , SURVL 

COMMON /C ALCON/ACTWT , ACTCMG, AC TLAM, MB IM1 ,MB I PI ,MB0M1 , MB0P1 , MMMl , 

1 HM1,HT,0TLR, DM LR, PITCH, CP, EXPON,TWW,CPTIP, TGROG , TB I , TBO , L AM8U A , 

2 TWL,ITMIN,ITMAX,NIP,IMS(2),BV(2),MV( 100 ), I V C 101), ITVI 100,2) , 

3 TV( 100,2) ,DTDMV (100,2 ),BETAV( 100,2) ,MH( 100, 2 I , DTDMH 1 100, 2 ) , 

4 BET AH ( 100,2) ,RMH( 100, 2) ,BFHl 100,2 ) ,RM( 100 ) , BE 1 1 00 ) , DBDM I 100 ) , 

5 SALI 100) , AAAI 100) 

COMMON /SLA/TSLI600) , U IN T( 6 ) 

C I MENS I ON MSL( 100) ,KK K ( 1 4 ) , P { 4 ) 

COMMON /SURVEL/ WTB ( 1 00 , 2 ) , WMB I 100 , 2 ) , XDOWN ( 400 ) , V ACROS ( 400 ) 

Cl MENS I ON W ( 2500 ) , BET A ( 2500) , DUDT ( 2500 ), DUDT T { 2500 ) , AAP(2500) , 

1 BBP12500) 

EQUIVALENCE (A, W), (All, 2), BETA), (A(l, 3) , DUDT ) , ( A ( 1 ,4 ) , DUDT T ) , 

1 (K» AAP ) , (RHO.BBP) 

INTEGER BLOAT, A ANDK, ER SOR , STRFN , SLCRD, SURVL , AATE MP , SURF , F I RS T , 

1 UPPER, SI, ST, SRW 

REAL K,KAK, LAMBDA, LMAX,MH,MLE* MR, MSL, MSP, MV.MVIMI 
CATA (KKK(J),J=4,14,2)/6*1H*/ 

C 

C CALL SLAVP AND SLAVBB THROUGHOUT THE REGION 
C 
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I T VU= I T V ( l » 1 ) 

I T VL = I TV ( 1,2) 

CO 10 I M= 1 , MB I M 1 
10 CALL SLAVP( IM f ITVU.ITVL) 

CO 20 I M= MB I , MBO 
I- 0 

20 CALL SL AVBB ( I M ) 

90 ITVU = I T V ( MBOP 1 * 1 ) 

ITVL = irV(MU0Pl«2) 

CO 100 I M=MBOP 1 1 MM 
100 CALL SLAVP(IM f [TVUfITVL) 

C 

C PLOT STREAMLINES 
C 

IF ( SLC RD . LE . 0 ) RETURN 
CQ 110 I M= 1 f MM 
110 M SL ( IM) * MV(IM) 

KKK ( 1) = 7 
KKK (2) = 6 
KKK ( 3 ) = MM 
P ( L ) = i. 

P ( 3 ) = 0. 

P ( 4 ) = 0. 

WR I rE(6, 1000) 

CALL PLOTMYI MSL,TSL, KKK, P) 

WR I TE ( 6 , 1010) 

RETURN 

1000 FORMAT ( 2HP T , 50X , 40HS TRE AML I NE PLOTS FOR REDUCED WEIGHT FLOW) 

1010 FORMAT (2HPL»40X « 7 OH STREAM LINES ARE PLOTTED WITH THETA ACRUSS THF 
1 PAGE AND M DOWN THE PAGE) 

END 


SUBROUTINE SLAV 
C 

C SLAV C^LCULATEj* RHO*W-SUB-M THROUGHOUT THE REGION AND ON THE BLADE 
C SURFACES, AND CALCULATES THE STREAMLINE LOCATIONS 
C 

COMMON SRW, ITER, IEND,LER(2) ,NER(2) 

COMMON / AUK RHO / A ( 2500 « 4 ) f U ( 2500 ) , K ( 2500 ) , RHO l 2 500 ) 

COMMON / INP/GAM, AR,T I PtRHOIP , WTFL , OMEGA , OR F , BE TA I , BE T AO , REDF AC , 

1 DENTOL ,MBI , MBO, MM , NBB I , NBL , NRSP , MR ( 50 ) , RMSP ( 50 ) , BESP ( 50 ) , 

2 BLDAT, AANOK, ERSOR, STRFNtSLCRD, INTVL,SURVL 

COMMON /CAL CON/ ACT WT , ACT CM G, AC TL AM » MB I M 1 , MB I P 1 ,M60M1 , MBO Pi , MMM1 , 

1 HMl ,HT,DTLR,DMLR, PITCH»CP, EXPON ,TWW»CPTIP*TGROG,TB I, TBO, LAMBDA, 

2 TWL , I TM IN , I TMAX ,N I P , I MS ( 2 ) , BV < 2 ) • MV ( 1 00 ) , l V ( 10 1 ) , I TV( 100 , 2 ) , 

3 TVC ICO, 2 ) , DTDMV ( 1 00 , 2 ) , BE TA V ( 100, 2 ) , MH( 100 , 2 ) , D TDMH ( l 00 , 2 ) , 

4 BET AH ( 100,2) , RMH( ICO, 2) , B EH { 100 , 2 ) , RM ( 1 00 ) , BE ( 100) , DBDM ( 100 > , 

5 SAL ( 100 ) , AAAI 100) 

COMMON /SLA/TSL (600) ,UINT( 6) 

C I MENS I ON TSP ( SO ) ,USP ( BO ) , T INI (6) ,DDT( 50 ) 

COMMON /SURVEL/ WTB ( 1 00 , 2 ) , WMB ( i 00 , 2 ) , XDOWN ( 400 ) , YACROS { 400 ) 

C I MENS ION W ( 2t> 00 ), BETA ( 2500 ) , DUDT ( 2500 ) , DUDTT { 2500 ), AAP ( 2500 ) , 

1 BBPI2500) 

EQUIVALENCE ( A , W ) , ( A ( 1 » 2 ) « BETA ) » ( A ( 1 » 3 ) , DUDT ) , ( A { 1 , 4) , DUDTT ) , 

1 (K, AAP ) , ( RHO, BBP) 

INTEGER BLOAT, AANDK,ERSOR, STRF N , SLC RD , SURVL , AATEMP , SURF , F I R ST , 
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1 UPPER, SI, ST, SRW 

REAL K,KAK, LAMBDA, LMAX,MH,MLE, MR, MSL, MSP, MV, MVIMl 

C SLAVP CALCULATES ALONG VERTICAL MESH LINES WHICH DO NOT 
C INTERSECT BLADES 
C 

ENTRY SLAVP ( IM, ITVU, ITVL) 

LOC= 0 

NSP= ITVL-ITVU+2 
IP = I V (I M ) - I 
CO 10 I T= l , NSP 
IP = IP+l 

TSP(IT) = FLOAT! IT+I TVU- 1 ) *HT 
10 LSP ( IT )= U( IP) 

LSP(NSP) = USP(1)+1. 

IP = IV(IM) 

INTU = INT (U ( I P ) *f>, ) 

IF { U ( IP) ,GT.O. ) I NTU = I NTU+ 1 
CO 20 J=1 , 5 

UINT(J) = FLOAT ( INTU ) /5 . 

20 INTU = INTU+l 

L I N T ( 6 ) = UINT(l) 

GO TO 100 
C 

C SLAVBB CALCULATES ALONG VERTICAL MESHLINES WHICH INTERSECT BLADES 

ENTRY SLAVBB ( I M ) 

LOC = 1 

ITVUP1 = ITV(IM»1) 

ITVLM1 = ITVIIM.2) 

ITVU = ITVUP1-1 
ITVL = I TVLM1+1 
ASP = ITVL-ITVU+1 
TSP(l) = TV ( I M , 1 ) 

TSP(NSP) * T V ( IM , 2 ) 

USP{ 1) = BV(1) 

LSP(NSP) = B V ( 2 ) 

IP = I V ( I M ) - 1 

NSPM1 = NSP-1 

IF ( 2.GT .NSPM1 ) GO TO 70 

CO 60 IT=2,NSPM1 

IP = IP+l 

TSP(IT) = FLOAT ( IT+ITVU-1)*HT 
60 USP( IT) = U( IP ) 

70 CO BO 1=1,6 

80 LINT(I) = FLOAT ( I-,l ) / 5 . 

c i 

C FOR BOTH SLAVP ANO SLAVBB, CALCULATE RH0+W-SU8-M IN THE REGION, AND 
C RHO*W AT VERTICAL MESH LINE INTERSECTIONS ON THE BLADE SURFACES 

C s 

100 CONTINUE I 

CALL SPLINE(TSP,USP, NSP,DDT,AAA) 

IT = LOC 
IPU = I V I I M ) 

IPL = I V ( I M+l )-l 
CO 110 I P= IPU, IPL 
IT = IT + 1 

CUDri IP) = DDT ( IT) 

110 CUDTTI IP) = AAAI IT) 

120 IF (LOC.EQ.O) GO TO 1 TO 
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ViMft I IM « 1 ) = DOT ( IMWTFL/BEI IM)/RM( IM) 

W MB I I M * 2 ) = DDTCNSP ) *WTFL/BE( IMJ/RMI IM) 

RMDTU2 = l RM ( IM) *DTDMV ( IM, 1 ) J **2 
RMU1L2 = I RM ( IM)»DTDMV( IM»2) )**2 
IF (RMDTU2.GT . LOOOO. ) WM8(IM,1) = 0. 

IF ( RM0TL2.GT. 10000. ) WMB(IM,2J = 0. 

WMBIIM, 1) = ABSIWMBl IM,1))*S0RT(1. +RMDTU2 ) 

WMBIIM, 2) = ABSIWMBl IM»2))*SQRT( I.+RMDTL2) 

130 IF ISICRD.LE.O) RETURN 
M = 6 

CALL SPLINT (USPtTSPtNSPfUINTtNl .TINT, A A A) 

CO 140 J= 1 » 6 
L = ( J- 1 ) * M M+ I M 

140 TSL I L ) = TINT I J) 

IF (IM.EQ.l) WRITEI6, 1000) 

WRITEI6, 1010) MV I IM) , (UINT(J)»TINTIJ).J=1»6) 

R E TURN 

1000 FORMAT I IH 1 30X , 46HSTR EAML INE COORDINATES FOR REDUCED WEIGHT FLOW/ 
1 1HL , 4X , 1 2HM COORDI NATL, 31 7X, 10HSTREAM FN . , 10X , 5HTHETA , 4X ) // ) 
1010 FORMAT! IX, 70 18. 7/1 19X.6G 18.7) ) 

END 


SUBROUTINE TANG 

TANG CALCULATES RHO*W-SUB-THETA AND THEN RHG*W THROUGHOUT THE REGION 
AND ON THE BLAuE SURFACES, AND CALCULATES THE VELUCITY ANGLE, BETA, 
THROUGHOUT THE REGION 


C 

C 

C 


COMMON SRW, I TER, I END , LER (2 ) , NERI 2 ) 

COMMON / AUKRHG/ A I 250C . 4 ) , U ( 2500 ) , K I 2500 ) , RHO! 2500 ) 

CCMN'ON /INP/GAM,AR»T IPfRHOIPfWTFLf OMEGA f ORF v BE T A I « BE TAG t REOF AC ? 

1 OENTCiL ,MBI ,MBO,MM , NB8 I, NBL, NRSP, MR ( 50 ) , RMSP ( 50) ,BESP( 50) , 

2 4LDAT, AANDK,ERSOR* STRFN,SLCRD, INTVL.SURVL 

COMMON /C ALCON/ AC TWT , ACTGMG, ACT LAM ,MB I M 1 , MB I P 1 ,MB0M1 , MBUP1 , MMM1 , 

1 HMi,HT,DTLR, DMLR, PITCH, CP, EXPON,TWW,CPTIP , TuROG, TBI ,T 80, LAMBDA, 

2 TWL , ITMIN, I TMAX , N I P , IMSI2),BV(2),MV(100),IV(10I),ITV(100,2), 

3 TV! 100,2) » DTDMV 1100,2 ) , BE TAV t 100,2) ,MH( 100,2) ,DTDMHl 100,2) , 

4 BET AH I 100 , 2 ) , RMH! ICO, 2) , BEH ( 100, 2 ) , RM 1 100 ) ,bE ( 100 ) , OB DM ( 100) , 

5 SALI 100) , AAAI 100) 

COMMON /SURVEL/ WTB I 1 GO , 2 ) , W MB I 100 , 2 ) , XDOWN 1 400 ) , Y ACROS I 400 ) 

COMMON /RHOS/RHOHB (100,2), RHCVB l 100,2) 

DIMENSION W ( 2500 ) , BET A ( 2 500 ) , DUDT I 2500 ) , DUDTT ( 2500 ) , AAP ( 2500 ) , 

1 B3PI2500) , „ 

EQUIVALENCE I A , W ) , ( A I 1 , 2 ) « BE TA ) , I A I 1 , 3 ) ,DUDT ),(A(I,4),DU0TT), 

1 (K, AAP ) , (RHO, BBP) 

C I MENS ION SPMI 100 ) ,USP( ICO ) , DDT I 100) ,DUDM( 100) ,DUDMM I 100 ) , 

1 UUDTMUOO) 

C I MENS I ON DWDMI 100) ,W IP! 100) CIDCT 

INTEGER BLOAT , A ANOK, ER SOR, STRFN , SLCRD , SURVL , A A TEMP , SURF , F I RST , 

1 UPPER, SI, ST, SRW 

REAL K » K AK » L AMBD A , LM AX , MH» ML E» MR » MSL « MSP ,MV» MV IM l 

EXTERNAL BL1.BL2 

PERFORM CALCULATIONS ALONG ONE HORIZONTAL LINE AT A TIME 
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IT = ITMIN 

10 IF (IT.GT.ITMAX) RETURN 
SI = 0 
C 

C ON THE GIVEN HORIZONTAL MESH LINE, FIND A FIRST POINT IN THE REGION 
C 

I F ( IT.GE.O.ANO. IT.LT.NBBI ) GO TO 60 
IM = MB I M 1 
20 I M= IM+1 

I F ( IM.GT.M8U) GO TO 200 
SURF = 1 

IF( IT.GE. I TV l IM,1) .AND. IT.LT.ITVI IM-1,1 ) } GO TU 70 
IFIIM. EQ. MROP1.ANO.IT. EQ.ITVIM BO »l J-l. AND. I TVIMBO, 1 ) - I TV ( MBO, 2 I 
l +NBBI.EQ.2) GO TO 70 
SURF = 2 

IF (I T. L E. I TV ( I M, 2 ) . AND. I T.GT . I TV I I M-l , 2 J ) GO TO 70 
GO TO 20 

FIRST POINT IS ON BOUNDARY A-H 

60 I M 1 = I 
IM = I 

SPM (1) = MV ( 1 ) 

LSP (1) = U( IT+l ) 

GO TO 90 

FIRST POINT IS ON A BLADE SURFACE 

70 SI - SURF 
I M 1 = I M— I 
I M2 = IM 

TH = FLOAT! IT)*HT 
M VI Ml = MV ( IMi ) 

IF ( IM . EQ . Mu IP 1 ) MVIM1 = MV I Ml + { MV ( I M2 I-MV IM I ) / 1000. 

LER ( 2) = 5 

BLCO (VIA ROOT) CALL NO. 5 

IF ( SI. EQ. 1 . AND. IM1.NE . MBO) CALL ROOT! MV IMI , MV ( I M2 I , TH, BL 1 , DTL R, 

1 ANS.AAA) 

LER12) = 6 

BLCD ! V I A ROOT) CALL NO. 6 

I F ( S 1. EQ. 2 ) C ALL ROOTfMVIMl , MV { I M2 ) , TH , BL2 , DTLR , ANS , AAA ) 

IF! SI.EQ. 1. AND. IMI. EQ. MBO) ANS = MV(MBG) 

SPM(IMl) = ANS 
USP ! IMI ) - BVtSl) 

MOVE ALONG HORIZONTAL MESH LINE UNTIL MESH LINE INTERSECTS BOUNDARY 

90 IF! IM.LT.MBI.OR. IM.GT.MBU) GO TO 120 
SURF = 1 

IF! IT.LT.ITVIIM, SURF) .AND. IT.GE. I TV ( IM-1, SURF) ) GO TO 1A0 
SURF = 2 

IF ( IT.GT. ITV! IM, SURF ) . AND. IT.LE. ITV! IM-l .SURF) ) GO TO 1A0 
120 S PM ( IM ) = MV(IM) 

IP = I PF ( IM, IT) 

LSP! IM) = UUP) 

IF UM.EQ.MM) GO TO 1 30 
I M= IM 1 
GC TO 90 
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FINAL POINT IS ON BOUNDARY D-E 

130 IMT = MM 
GO TO 150 

FINAL POINT IS ON A BLADE SURFACE 

140 ST ^ SURF 
I M T = I M 

IMT«1= IMT-l 
TH = FLOAT ( IT) *HT 
MVIMl = I MV( IMTMl ) 

IF ( IM.EQ.MBIPl) MVIMl = M V I M 1 + I MV ( I M T ) -M V I M I ) / 1000. 

LER12) = 7 

ELCI) (VIA ROOT) CALL NO. 7 

IFUST.EQ. l.AND. IMT.NE.MBDCALL ROOT (MVIMl » MV ( I M T ) , TH , BL l * 

l DTLR * ANS » AAA ) 

LERI?) = 8 

BLCO (VIA ROOT) CALL NO. 8 

IF ( ST. EC. 2 ) CALL ROOTtMVIMl , MV ( I MT ) , TH , RL 2 , OTLR , ANS . AAA ) 

IF (ST.EO. 1 . AND. IMT .EQ.MB I ) ANS = MV(MBI) 

SPM(IMT) = ANS 
USP ( IMT ) = RV(ST) 

150 NSP= IMT-IMl+l 

CALL SPLINE (SPM(IMl), USP ( I M l ) , NSP , DUDM ( I M l ) . DUDMM I I M 1 ) ) 

CALCULATE RH0*W ON THE BLADE SURFACES 

F I R ST= 1 
L AS 1= MM 

I F ( IM1.FQ.1) GO TO 160 
FIRST = I M2 

CALL SEARCH (SPM( I Ml ) »S1»IHS) 

ANS =-DUDM( I M 1 )*WTFL/BEH( IHS.S1 ) 

WTB ( IHS » S 1 ) = ABS( ANS )*SCRH 1. + 1. /(RMH( IHS.Sl ) *OTDMH( IHS , SI ) ) **2 ) 

CDT(IMl) -=-DUDM ( IMl)/CTDMH( I HS , S 1 ) 

V» I P ( I M 1 ) = nTb ( I HS » S 1 ) / RHOHBl IHS.Sl ) 

160 IF( IMT.CQ.MM) GO TO 170 
LAST = IMTMl 

CALL SEARCH ( SPM ( I MT ) , ST , I HS ) 

ANS =- DUD M ( IMT)*WTFL/BEH(IHS,ST) 

WTB( IHS»ST) = A6 S ( AN S ) »SGR T ( 1 . + 1 . / ( RMH ( IHS.ST)*DTDMH(IHS.ST) )«»2) 
C D T ( IMT) =-DUDM( I MT) / CTDMH ( I HS , ST ) 

K I P ( IMT) = wTBUHS.ST) / RHOHB ( IHS * ST ) 

CALCULATE RHO*w-SUB-THETA AND THEN RHO.W AND BETA IN THE REGION 

170 IF (FIRST. GT. LAST ) GO TO 190 
CO L 80 I = FIRST,LAST 
IP = I P F ( I, IT) 

CDT ( ( ) = DIJOT (IP) 

RWM = DDT ( I ) /RM( I ) 

RWT =-DUDM ( I ) 

MIP) = SORT ( RWT**2+RWM»*2 )/BE( I )*WTFL 
TWLMR = 2 . «-UMEGA *L AM BE A- (OMEGA*RM ( I ))»*2 
LER(l) = 5 
C CENSTY CALL NO. 5 

CALL DENSTYI W ( I P ) . RHO ( IP ) » ANS , TWl.MR , CPT I P » EXPUN > RHO IP,GAM,AR,TIP) 
V% t I ° ) = ANS 
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WIPI I) = W( IP) 

BETMIP) = A T AN2 I RWT ,RWM)*57.295779 
180 CONTINUE 

IF( IEND.LT.O) GO TO 1<30 

CALL SPLINE I SPM ( I Ml ) , CD T( I M 1 ) , NSP ,DUDTM ( I Ml ) , AAA ( I M I ) ) 

CALL SPLINE I SPM ( IM1 ) , W I P( I M 1 ) , NSP ,DWDM ( IM 1 ) , AAA I I M 1 ) ) 

DO 185 I=FIRST,LAST 
IP = I P F ( I , IT) 

SBETA = SINIBETAI IP)/57. 295779) 

CBt TA = SQRTI l.-SBETA**2 ) 

A A P ( IP) = SbETA**2*( 2 .*DUDTM( I )/DUDM( I )-OUDT ( IP ) /OUDMI I ) **2* 

1 UUDMMf I )-DUDTT( IP )/CUDT( IP) ) + S AL I I )*SBETA/CBETA*( l.+C8ETA**2) 
BBP ( IP ) = RM ( I ) /CBET A* I 2 .*ACTOMG*SAL 1 1 ) + SBETA* DWDM I I ) /RED FAC ) 

185 CONTINUE 
190 CONTINUE 

I F { IMT.NE.MM) GO TO 20 
200 IT = IT+I 
GO TO 10 
END 


SUBROUTINE SEARCH ID I ST , SURF , [ S > 

SEARCH LOCATES THE POSITION OF A GIVEN VALUE OF M IN THE MH ARRAY 


COMMON /CALC ON/ ACTWT , ACT CMG, ACTLAM , MB I Ml ,MB I PI , MBUM 1 , MB0P1 , MMM1 , 

1 HMl,HT,DTLR,DMLR, PITCH, CP, CXPON,TWW,C PT IP, TGROG.TB I ,TBO, LAMBDA, 

2 TWL , ITMINt ITMAX.NIP, IMS ( 2 ) » BVI 2) ,MV( LOO) , IV ( 101 ) , ITVI 100*2) , 

3 TV! 100,2) ,DTDMV I 100,2 ), BETA VI 100,2) ,MH( 100,2) .DTDMHI 100,2) , 

A BET AH ( 100,2 ) ,RMH{ 1 CO, 2 ) , B EH ( 100 , 2 ) , RM { 100 ) , BE ( 1 00 ) , DBDM ( 100 ) , 

5 SAL ( 100 ) , AAAI 100) 

INTEGER flLDAT, AANDK, ER SO R,STRFN,SLCRD,SURVL,AA TEMP, SURF, FIRST, 

1 UPPER, SI, ST, SRW 

REAL K, KAK, LAMBDA, LM AX, MH,MLE,MR,MSL, MSP, MV, MVIMI 
CO 10 [=1,100 

IF l ABS (MH ( I , SURF ) -D I ST ) .GT. DMLR ) GO TO 10 
IS = I 
RETURN 
10 CONTINUE 

WRITE (6,1000) DI ST, SURF 
STOP 

1000 FORMAT I 38HL SEARCH CANNCT FIND M IN THE MH ARRAY/ 7H DIST =,G1A.6, 
110X.6HSURF = , G l A . 6 ) 

END 
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SUBROUTINE VELOCY 


VELCCY CALLS SUBROUTINES TO CALCULATE DENSITIES AND VELOCITIES 
THROUGHOUT THE REGION AND ON THE BLADE SURFACES, AND IT PLOTS 
THE SURFACE VELOCITIES 

COMMON /AUKRHC/ A ( 2500 , A ) , U ( 2500 ) , K ( 2500 ) , RHOI 2500 ) 

COMMON /INP/GAM, AR, T I P , R HO IP , WTFL , OMEGA , ORF , BET A I , BE T AO, RE OF AC , 

1 DENTOL , MBI , M80.MM , NRB I , N8L , NRSP ,MR ( 50 ) , RMSP ( 50 ) , BESP C 50 ) , 

2 HLDAT, AANDK » ERSOR » STRFN.SLCRD, INTVL,SURVL 
COMMON /CALCON/ ACTWT , ACTCMG , AC TL AM , MB I M 1 , MB I P 1 , MB0M1 , MBOP 1 , MMM l , 

1 HMI ,HT , DTLR ,DMLR, P I TC H, C P , E XPON , TWW , CPT I P , TGROG , FBI »T80, LAMBDA, 

2 TWL. ITMIN, ITMAX.NIP, IMS( 2 ) ,BV(2 > ,MV( 100) , IV ( 101 ) , ITVl 100,2) , 

3 rv( 100,2) ,UTDMV(100,2),BETAV( 1 00, 2 ) , MH ( 1 00 , 2 ) , DTDMH (100,2), 

A BET AH ( 100,2) , RMH( 1 CO, 2) , BEH t 100, 2 ) , RM 11 00 ) , BE I 1 00 ) ,DBDM ( 1 00 ) , 

5 SAL ( 100) , AAAI 100) 

DIMENSION KKK ( LA) 

COMMON /SURVEL/ WTB ( l 00 , 2 ) , W MB I 100 , 2 ) , XDOWN ( AOO ) , YACROS ( AOO ) 
DIMENSION W( 2500) .BETA (2500) ,DUDT{ 2500) .DUDTTI 2500) ,AAP( 2500) , 

1 BBPI2500) 

EQUl VALENCE ( A, w ) , ( A ( 1 , 2 ) , BETA ) , ( A ( 1,3) ,DUOT) , I A ( 1 , A ) , DUDTT ) , 

1 (K,AAP) , (RHO.BBP) 

INTEGER BL DAT, AANDK, E RSOR, S TRFN , SLCRD , SURVL , AATEMP , SURF , F I RST , 

1 UPPER, SI, ST, SRW 

REAL K,KAK, LAMBDA, LMAX.Mh, ML F,MR,MSL,MSP, MV, MVIM1 
CAT A KKKIAl/LH*/, KKK I 6 ) / 1H0/ , KKK I 8 ) / LH+ / ,KKK ( 1 0 ) / IHX/ 

CALL VFLP , VELBB, AND VELSUR THROUGHOUT THE REGION 

I F { [ NT VL . GT . 0 ) C ALL V EL P ( 1 , MB I M l ) 

CALL VELRBIMBI ,MBO) 

20 IF! INTVL.GT.01CALL VELP I MBOP 1 , MM ) 

CALL VELSUR 

PREPARE INPUT ARRAYS FOR PLOT OF VELOCITIES 

I F ( SURVL. LE. 0 ) RETURN 
NP2 = 0 
C TANGENTIAL COMPONENTS 
CO 50 S UR F = 1 , 2 
NP1 = NP2 
IMSS = IMS ( SURF ) 

IF ( IMSS.LT. i ) GO TO AO 
CO 30 IHS= 1 » IMSS 

IF (ABSIDTDMHI I HS , SUR F ) *RMH ( IHS, SURF) ).LT.. 57735) GO TO 30 
NPI = NPl+I 

YACROS(NPl) = WTB( IHS, SURF ) 

XDOwN ( NPI) = MH I IHS, SURF ) 

30 CONTINUE 

AO KKK ( 2»SURF+1 ) = NP1-NP2 
50 NP2 = NPI 

C MERIDIONAL COMPONENTS 
CO 80 SURF= I » 2 
NPI = NP2 

CO GO I M=MB I PI , MB0M1 

IF (ABSIDTDMVI IM, SURF )*RM(IM) ) .GT. 1.7321) GO TO 60 
NPI = NPl+l 

YACROS (NPI) = WMBIIM.SURF) 

XDC)wN( NPI ) = MV( IM) 

60 CONTINUE 
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70 KKK { 2*SURF+5 ) = NP1-NP2 
80 KP2 = NP1 
C 

C PLOT VELOCITIES 
C 

KKK ( 1 ) = 1 
KKK ( 2 ) = A 
P = 5. 

WR I TE ( 6 , 1000) 

CALL PLOTMY(XDOWN,YACROS,KKK,P) 

WRITEC6.101G) 

RETURN 

1000 FCR V AT( 2HPT, SOX, A8HBI APE SURFACE VELOCITIES FOR REDUCED WEIGHT FLO 
1W ) 

1010 FORMAT (2HPL,37X,63HVELOCITY(W) VS. MERIDIONAL STREAMLINE DISTANCE 
l(M) DOWN THE PAGE / 2HPL / 

2 2HPL » 5 OX » 50H + - BLADE SURFACE 1, BASED ON MERIDIONAL COMPONENT/ 

3 2HPL , 50X , 50H* - BLADE SURFACE 1, BASED ON TANGENTIAL COMPONENT/ 

A 2HPL » 50X . 50HX - BLAOE SURFACE 2, BASED ON MERIDIONAL COMPONENT/ 

5 2HPL , 50X , 50H0 - BLADE SURFACE 2. BASFD ON TANGENTIAL COMPONFNT) 

END 


SUBROUTINE VEL 
C 

C VEL C Al.CUL AT ES DENSITIES AND VELOCITIES FROM THE PRODUCT OF 
C DENSITY TIMES VELOCITY 
C 

common SR W, I TER. I END ,LER (2 ) , NER{ 2) 

COMMON / AUKRHO/ A { 2500 , A ) , U ( 2500 ) , K ( 2500 ) , RHO I 2500 ) 

COMMON /INP/GAM, Art, T IP. RHO I P , WTFL , OMEGA , ORF , BETA I ,F,ETAO, REDFAC , 

1 DENT 0L , MB I , MOO, MM. NBBI.NBL ,MRSP f MR( 50) , RM SP I 50 ) , BFSP ( 50 ) , 

2 BLOAT . AANDK. ERSOR . STRFN » $ LCRD » I N TVL , SUR VL 
COMMON /CALCON/ ACTWT , AC TCMG , ACTl AM, MB I Ml , MB I P 1 , MB0M1 , MBOP 1 , MMM l , 

1 HM1 .HT.DTLR.DMLR, P I TCH , CP , EXPON , TWW , CPT I P , TGROG , TO I , TBO, LAMBDA, 

2 TWL, I TM IN, I TMAX.NI P, IMS < 2 ) , BV( 2 ) ,MV( 100) , IV ( 101) , I T V C LOO, 2 ) , 

3 I V( 100,2) , OTDMV { 100,2 ), BETA V( 100,2) ,MH{ 100,2) »DTDMH ( 100,2 ) , 

A BET AH ( 100,2) ,RMH{ 100,2) ,BEH( 100,21 ,RM( 1 00 ),BE(100),DBDM( 100), 

5 SALI 100) , AAAI 100) 

COMMON /RH0S/RH0H8 ( 100,2 I.RHOVBI 100,2) 

Cl McNS ION WWCKMI 100, 2 ) , WrtCRT I 100,2) , SUREL ( 100,2) 

COMMON /SURVEL/ WTBI 100,2) » WMB { 100,2) .XDOWN(AOO) .YACROSI AOO) 

Cl Mf NS ION W( 2500) .BETA (2500) , DUUT ( 2500 ), DUDTT ( 2500 ) , AAPI 2500) , 

1 EBP ( 2500 ) 

EQUIVALENCE (A, W), (All, 2), BETA), (All, 3) , DUDT ), (All, A), DUDTT), 

1 (K, AAP ) , ( RHO, 8BP) 

VELP CMCULATES ALONG VERTICAL MESH LINES WHICH DO NOT 
INTERSECT BLADES 

ENTRY VELP ( F IRST , LAST ) 

INTEGER BLOAT, AANDK, ERSOR, STRFN, SLCRD , SUR VL , AATEMP , SURF , F I RST , 

1 UPPER, SI, ST, SRW 

REAL K,KAK, LAMBDA, LMAX, ME, MLF,MR»MSL»MSP»MV»MvlMl 
IF ( F IRST. GT. LAS T ) RETURN 
IF (FIRST. Ew.l) WRITE (6,1000) 

CO 20 IM = F IRST, LAST 
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IPU = IV( I M ) 

I PL = IPU + NdBI-i 

WRITE (6,1010) IM, (W( IP) ,BETA( IP) , IP=IPU,IPL> 

20 CONTINUE 
RETURN 

VELBR CALCULATES ALONG VERTICAL MtSH LINES WHICH INTERSECT BLADES 


C 


c 


c 


c 


c 


ENTRY VELGR(FIRST,LAST) 

IF ( FIRST. GT. LAST ) RETURN 
IF (FIRST. NE. MBI ) GO TO 30 
RELER = .0 
SURFLI MBI , 1 ) = 0. 

SURFL ( MB I , 2 ) = 0. 

30 CO JO I M=F IRST , L AST 
I TVU = I TV l IM, 1) 

ITVL = I TV ( I M , 2 ) 

I PUP 1 = I PF ( I M » I TVU) 

IPLM1 = I PF ( IM, ITVL) 

TWLMR = 2 . *GMEG A *L AM B 0 A- (OME G A*RM ( I M ) ) * * 2 
WCR = SQRT ( TGROG*T IP* 1 1 .-TWLMR/CPT IP ) ) 

IF ( ITVL.LT. ITVU) GO TO 50 
ALONG THE LINE BETWEEN BLADES 
IF ( INTVL.LE.O) GO TO 50 

WRITE (6,1010) IM, (W( IP) .BETA! IP) , IP=IPUP1. IPLM1) 


ON THE UPPER SURFACE 
50 PHOB = RHOVB ( IM, 1 ) 

LER( 1)=6 

CALL DENSTY ( WMB ( IM,1 ) ,RHCVB( IM,1) , AN S , TWLMR ,CP T l P , FXPON, RHO IP, 
1 GAM, AR, TIP) 

WMB ( IM , l ) = ANS 
WWCRM(IM,1) = w'MBt IM, 1 )/WCR 
IF ( IM. EQ. MB I ) GO TO 6C 
CELTV = TV( IM-1, 1 )-TV ( IM.l ) 

SURFL ( I M , 1 ) = SURFLI IM-1, 1J+SQRTI (MVI IM)-MV(IM-l) )**2 + 

1 (DELTV*(RM( IMI+RMI IM-1) )/2. )**2) 

60 RELER = AMAx 1 ( RELER, ABS ( (RHOR-RHOVB l I M, 1 )) /RHOVB ( I M , l )) ) 

ON THE LOWER SURFACE 
RH03 = RHOVB ( l M, 2 ) 


L ER ( 1 ) =7 

CENSTY CALL NO. 7 

CALL DENSTY(WMB( IM,2) , RHOVB ( IM,2),ANS»TWLMR,CPTIP»EXP0N,RH01P» 
1 G AM , AR , T I P ) 

WMB ( I M , 2 ) = ANS 

WWCRMI I M, 2 ) = WMBI IM.2J/WCR 
I F ( IM. EQ. MB l ) GL) TO 70 
CELTV = T V ( I M- 1 , 2 )-T V ( I M , 2 ) 

SURFLI I M , 2 J = SURFL( IM-1,2 ) + SQRT( ( MV ( IMJ-MVI IM-L) )**2 + 

1 (OELTV* I RM( IM ) +RMI IM-1) )/ 2. ) **2 ) 

70 RELER = AMAX1 ( RELER, ARS( (RHOR-RHOVB ( IM»2) ) /RHOVB ( I M , 2 ) ) ) 

RETURN 


C VELSUR CALCULATES ALONG A BLADE SURFACE 
C 

ENTRY VELSUR 
CO 70 SURF =1,2 
I MSS = IMS(SURF) 

I F ( IMSS.EQ.O) GO TO 9 0 
CO 30 I H S = 1 » I MSS 
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TWLf*R = 2.*0MEGA*LAMBCA- (OMEGA*RMH| IHS, SURF) )**2 
UCR = SQRT(TGROG*TIP* { 1 . -TWLMR /CP T I P ) ) 

RHOH = RHOHB (IHS, SUR F ) 

LER ( 1 ) =8 

C CENSTY CALL NO. 8 

CALL DENSTY ( WT8! IHS, SURF ) , RHOH3 l I HS , SURF ) , ANS , TWLMR ,CPT I P , 

1 FXPON, RHGIP,GAM,AR,TI P) 

WTB( IHS, SURF) = ANS 

VWCRT! IHS, SURF ) = WTB { IHS, SURF ) /WCR 
80 RELFR = AMAX1(RELER,ABS! (RHOfi-RHOHB ( IHS , SURF ) ) /RHOHB ( IHS, SURF ) H 
90 CONTINUE 

IF(RELER.LT.OtNTOL) I END = IEND+1 
WR I TE ( 6 , l 080 ) ITER.RELER 

WRITE ALL BLADE SURFACE VELOCITIES 

IF (SURVL.LC.O) RETURN 
WRITE(6,1020) 

WRITE (6, 10 AO) !MV( IM ) ,WMB( IM , I ) , BETAV ( IM, L ) , SURFL { IM, 1), 

1 WWCRMI IM, l ) ,WMB( IM,2) , BETAV ( 1 M , 2 I » SURFL ( IM, 2) ,wWCRM( IM,2) , 

2 IM=M8 I , Mt)0 ) 

WRITEI6, 1050) 

CO 100 SURF= l , 2 
I MS S = IMS t SURF ) 

IFIIMSS.LT. i) GO TO 100 
WRITE(6,1 060 ) SURF 

WR I TE ( 6 » l 070 ) ( MH (IHS, SURF ), WTB I IHS, SURF ), BETAH ( IHS , SURF ) , 

1 WWCRTI IHS, SURF) , IhS=l,IMSS) 

100 CONTINUE 
RETURN 

1000 FORMAT! IH1/40X.34HVEL0CI TIES AT INTERIOR MESH P0INTS/45X, 

1 23HF0R REDUCED WEIGHT FLOW) 

1010 FORMAT ( IHL,3H1M=, 13, 5 ( 24H VELOCITY ANGLF (DEG ) ) / 

1 !5X,5(G15.4»F9.2) ) ) 

1020 FORMAT { 1H1/I6X, 1H*,I8X,71HSURFACE VELOCITIES BASED ON MERIDIONAL C 
1CMP0NENTS - REDUCED WEIGHT FLOW , 1 8X, 1H*/ 16X, IH* , 53 X , 1H*5 3X , IH*/ 

2 16X, IH*, 19X, 15HBLACE SURFACE 1, 19X, 1H*,20X, 15HBLADE SURFACE 2, 

3 I8X, 1H*//X, 1HM.8X, IH* ,2 ! 3 X , 8HVELUC I TY , 3X , 23HANGLF ( DEG ) SURF. LE 
4NGTH,5X,5HW/WCR,6X,1H*) ) 

1040 FORMATULH ,G13.4,3H * , 2! G 1 2 . 4 , F9 . 2 , 2G 15. 4 , 3H *))) 

1050 FORMAT! IH 1/ 3X , 49HSURF ACE VELOCITIES BASED ON TANGENTIAL COMPONENTS 
I / 18X , 1 9 HR EDUCED WEIGHT FLOW) 

1060 FORMAT! //22X, I5HBLADE SURFACE , 1 1 / 7X , 1 HM, 10X , 8HVEL0C I T Y , 3X , 1 OHANG 
1LE t DEG ) , 3X,5Hw/WCR) 

1070 FORMAT I IH , 2G1 3 . 4 , F9 . 2 , G 15 . 4 ) 

1080 FORMAT ( 14HL ITERATION NO . , I 3 , 3X , 36HMAX I MUM RELATIVE CHANGE IN DENSI 
1TY =,G1 1.4) 

END 
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SUBROUTINE IVELCY 

COMMON SRW.ITfcR, IEND,LER(2) , NER(2) 

COMMON / AUKRHO/ A ( 2500 , 4 ) , U ( 2500 ) , K ( 2500 ) , RHO( 2500 ) 

COM WON / I NP/GAM » AR*T l P • RHO I P » W TFL * OMEGA t GRF * 8E TA I * BE TAGt REDFAC » 

1 QENTOL.MBI , MBC , MM , NBB I , NBL , NRSP , MR ( 50 ) , RMSP ( 50 ) ,0ESP(50) , 

2 BLOAT, AANOK.ERSOR, STRHN.SLCRO, INTVL,SUKVL 

COMMON /CALCON/ACTWT , ACT0MG,ACTLAM»MBIM1,MBIP1 ,MBOMl ,MB0P1 , MMMl » 

1 HMl,HT,DTLR,DMLR, PITCH, CP , EXPON » T MW, CPT I P , T GROG, TB I , T BO, LAMBDA , 

2 TWL, ITMIN, ITMAX,NIP, I MS ( 2 ),BVI2 ) ,MV( 100 ) , IV( 101 ) , l TV ( 100 , 2 ) , 

3 TV( 100 , 2 ) ,0 TOM VI 100,2 ), BETAVl 100,2 ) ,MH( 100, 2 ) ,l)TOMH( 100,2 ) , 

4 BET AH ( 100,2) ,RMH( 100, 2) , BEHl 100,2) ,RM( 100 ) , BE ( 100) ,DBDM( 100) , 

5 SAL( 100) , AAAI 100) 

COMMON /WWCRM/WWCRMl ICO, 2) , LABEL ( 100) 

COMMON /SURVEL/ WT B( 100, 2) , WMB I 100, 2 ) , XOOWN 1400) »YACR0S(400) 

C l MENS I ON W ( 2500 ), BETA ( 2 500 ) ,DU0T { 2500) , DUDTTI 2500) ,AAP(2500) , 
l B BP I 2500 ) 

EQUIVALENCE ( A,W ) , ( AC 1,2), BETA), (All, 3) , OUOT ) , ( A( 1 , 4 ) , DUOTT ) , 

1 ( K, AAP) , (RHO.BBP) 

C I MENS I ON KKK ( 14 ) 

INTEGER BLOAT, AANDK, ERSOR,STRFN, SLCRD.SURVL, AATEMP, SURF, SURFBV, 
1FIRST, UPPER, UPPRBV, SI ,ST ,SRW 
REAL K.KAK, LAMBDA, LMAX,MH,MLE,MR,MSL. MSP, MV, MVIM1 

LAMBDA = ACTLAM 

IF ( 1NTVL.GT.0 ) W R I TE ( 6 , 1 COO ) 

IF(l.GT.MBIMl) GO TO 20 
CO 10 I M= l , MB I Ml 
10 CALL VELGRAl IM) 

20 CO 30 I M= MR I P 1 , M30M1 
30 CALL VELGRR ( IM ) 

IF( MBOPl.GT.MM) GO TO 50 
CO 40 I M= MBOP I , MM 
40 CALL VE LGR A ( IM) 

IF(SURVL.LC.O) RETURN 

WRI TE16*1020) (MV ( IM) » WMB ( I M« 1 ) , WWCRM ( IM, L) , LABEL ( IM) ,WMR( IM,2) , 

1 WWCRM(IM,2)«LABEL( IM ) , I M = MB IP1 » MBOMl ) 

C PREPARE ARRAYS EUR PLOT CF VELOCITIES 
CO 60 IM = MB I P I , MBOM 1 
I = IM - MR I 
12 = I + MBOM 1 - MBI 
XDOWNI 1 ) = MV( IM ) 

YACROSl I ) = WMB ( I M , 1 ) 

60 YACROSl 12) = WMb(IM,2) 

KKK ( l ) = 0 

KKK ( 2 ) = 2 

KKK ( 3 ) = MBOM 1 - MBI 

P = l, 

C PLCT VELOCITIES 
WRITE (6,1030) 

CALL PLOT MY ( XDOW.N,YACROS, KKK , P ) 

WRITE (6,1040) 

RETURN 

1000 FORMAT ( IH 1 /40X » 34HVE LCC I T I ES AT INTERIOR MESH P0I.NTS/44X, 

1 27H ( BASED ON FULL WEIGHT FLOW)) 

1010 FORMAT ( 1H 1/ 16X , 1H* , 1 3X , 6 8HSURF ACE VELOCITIES BASED ON MERIDIONAL C 
I CM PONE NTS - FULL WFIGhT FLOW , 30X , 1H* / 16X , 1H* , 5 5X , I H* , 55X , 1H* / 1 6X , 
21H* ,20X , 15HBLADE SURFACE I , 20X , 1H*,20X, 1 5HBL ADc SURFACE 2, 20X,IH*/ 
37X , 1HM, 8X , IH*» 2 ( 3X» 8HVEL0C I T Y, 6X» 5HW/WCR , 33X , 1H* ) ) 

1020 FORMAT ( ( 1X,GL3.4, 3H » , 2 ( 2G 1 3 . 4 , A9 , 20X, 1 H* ) ) ) 
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1030 l F FULL T WFIGHr’FLUWn HBLADE SURFACE VFl -0C ITI ES/2HPT , 49X , 27H ( BA SED ON 

1040 FORMAT ( 2HP L , 3 7X , 6 3H V fcLOC I T Y ( W ) VS. MERIDIONAL STREAMLINE DISTANCE 
1(M) DOWN THE PAGE/2HPL/2HPL , 50X, 19H* - BLADE SURFACE l /2HPL . 50X . 1 9 
2H+ - BLADE SURFACE 2) 

END 


SUBROUTINE VELGRA(IM) 

COMMON SRW,IT£R,IEND,LER(2),NER(2) 

COMMON /AUKRHO/ A ( 2500 , 4 ) , U { 2500 ) , K ( 2500 ) , RHO C 2 500 ) 

COMMON / I NP/GAM , AR , T I P ,RFO IP , WTFL, OMEGA ( ORF, BE TAI , BETAO, RED FAC , 

1 OENTOL * Mb I t^B0fMM # NBBItNBLtNRSPfMR(50)tRMSP(50) * BE $P ( 50 ) 

2 BLDAT, AANDK,ERSOR, STR FN , SLCRD, I NTVL , SUR VL 
COMMON /CALCON/ACTWT, ACTCMG, ACTLAM,MBIMl f MBIPl,MBOMl ,MB0P1 ,MMM1, 


HMl »HTfDTLR»DMLR» P ITCH, CP ,E XPON , T WW > C PT I P , TGROG , TB I , T BO, LAMBDA, 
TWL, ITMIN, ITMAX.NIP, IMS(2),BVI2 J ,MV( 100), IV( 101 ) , I TV ( 100,2) , 

TV( 100, 2 ) ,l)TDMV (100,2 ) , BETAVI 100, 2) , MH( 100, 2) ,DTOMH( 100,2 ) * 

BET AH ( 100,2 ) ,RMH( 1 00, 2) , B FH ( 1 00 , 2 ) ,RM ( 100 ) , faE ( 1 00 ) , DBDM ( 1 00 1 , 
SAL( 100) , AAAI 100) * 

COMMON /D2T0M2/ D2TDM 2 ( 1 CO , 2 ) 

COMMON /WWCRM/WWCRMI 1 00 , 2 ) , L ABEL ( 100 ) 

DIMENSION WGRADI 50 ), THET A{ 50 ) , RWCB ( 50 ) , CBETA { 50 ) , A2 ( 50 ) , B2 ( 50) 
COMMON /SURVEL/ WT B{ 1 00 , 2) , WMR ( 100, 2 ) , XDOWN ( 4 IK) ) , YACROS { 400 ) 
DIMENSION W(2500) , BE T A ( 2 500 ) ,DUOT 12500) ,DUDTT(2500) ,AAP(2500) , 

1 EBP ( 2500 ) ’ 

EQUIVALENCE ( A ,W ) , ( A ( 1,2), BETA), (All, 3), DUDT ) , { A ( 1 , 4 ) , DUDTT ) , 

1 (K, AAP) , (RHO.BBP) 

INTEGER BLDAT, AANDK, E RSOR, STRFN, SLCRD, SURVL , AATEMP , SURF, SURFBV, 
1FIRST, UPPER, UPPRBV, SI ,ST,SRW 
INTEGER CHOKED 

REAL K»KAK, LAMBDA, LMAX,MH,MLF, MR, MSL, MSP, MV, MVIM1 
CATA CHOKED/6HCH0KED/ 


AORB = 1. 

IP = I V ( IM) 

WGRAD(l) = W ( I P ) /REDFAC 
IT1 = I TV ( IM, 1 ) 

IT = I T 1 
NSP = NBBI+1 
CO TO 10 


ENTRY VELGRbI IM) 

IP = I V ( IM)— 1 

WGRADI 1 ) = WM6 ( IM, 1 ) /REDEAC 
NSP = IV( I M+ 1 ) - I V ( IMJ+2 
AORB = 2. 

I T 1 = ITV( IM, 1 ) 

IT - IT 1-1 
10 NSPM 1 = NSP-l 

TORSAL = 2.*ACT0MG*RM( IM)*SALl IM) 

TWLMR = 2. *ACTOMG*LAMBCA-( ACTOMG*RM( IM) )«*2 
WCR = SORT ( TGROG*T IP* ( 1 . -TwLMR/CPT IP ) ) 
CELMAX = WCR/ 10. 

TOLERC = WCR/ 100. 

CO 20 1=1, NSPM1 

CBETA(I) = COS(BETA( IPI/57. 295779) 

THE TA ( I ) = HT*FLOAT( I T ) 
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A21 I ) = A AP ( IP) 

B2( I ) = BBPI IP) 

IT = IT+l 
20 IP = I P+1 

CBETA(NSP) = CBETA(l) 

A2INSP) = A2 ( 1 ) 

82 (NSP) = H2 ( 1 ) 

THETA(NSP) - HT *FL0AT (IT) 

IFIADRB.LE. 1. ) GO TO 30 

CBETA(l) = L./SQRT(1.+(RM( IM)*DTOMV( IM» 1 J )**2) 

SBETA1 = SORT ( l ,-CBFTA ( 1 ) **2 ) *S IGN l 1 . , OTDMV ( IM, 1 ) ) 

42(1) = (RM(IM)*CBETA(1) )**2*D2TDM2 ( IM , l ) + SAL ( I M ) * S3E TA1 / 

1 C8ETA( l )*( 1 . +CBET A ( 1) »*2) 

R2(L) = B2(2)+TCRSAL*( l./CBETAl 1)-1./C6ETA(2) ) 

THE T A( I) = TV( I M » l ) 

CBETA(NSP) = l./SQRT( 1 . ♦ (RM { IM ) *DTOMV ( I M ,2 ) ) **2 ) 

SBETAN * SORT ( l. -CBETAt NSP ) »*2 ) *SIGN( 1 . , DTOMV ( IM , 2 ) ) 

A2(NSP) = ( KM( IM ) *CBF TA ( NSP )) **2*D2TDM2 ( IM, 2 )+SAL ( IM )• SBETAN/ 

1 CBETA(NSP)*( l.+CBETAINSP) **2) 

P2 ( NSP ) = B2 (NSPM1 )+TORSAL*( l./CBETAINSP ) - 1 . /CBETA ( NSPM1 ) ) 
THETA ( N SP ) = TV ( IM,2) 

30 IND = 1 
AO CONTINUE 

DO 50 1=2, NSP 

WAS = WGRADt 1-1 )♦( A2 ( I-I )*WGRAO( 1-1 ) +B21 1-1 ) )* (THETA( I )- 

1 THETAd-l)) ... 

V>ASS = WGRAD ( 1-1 )♦( A2 ( I ) »WAS + B2( I ) >* (THETA ( I )-THETA( 1-1 ) ) 

50 WGRADt l ) = ( WAS+WASS ) /2. 

CO 60 1=1, NSP 

TTIP = l.-(WGRAD( I)** 2+T WLMR ) /CPT IP 
IFITTIP.GE. .0) GO TO 55 
WRITE (6,1010) IM 
WGRADt 1) = 0. 

WGRADt NSP ) = 0. 

GO TO 80 

55 RHOT = RHOIP*TTI P*»EXPON 
60 RWCb(I) = RHOT*WGRAD( I )*CBETA( I ) 

CALL INTGRL ( THETA, RWCB, NSP , AAA ) 

WTFI ES = BF ( IM ) *RM( I M ) »AAA ( NSP ) 

IF ( ABS( ACTWT-WTFLES ) .LE . ACT WT/ 100000. ) GO TO 70 

CALL CONTIN ( WGR AO ( 1 ) , WT FL E S , I ND, I M , ACTWT , DELMAX , TOLFRC ) 

IF ( IND.LT.6) GO TO AO 
IF ( IND.EQ.6) GO TO 65 
WRITE (6,1020) IM 
IF (AORB.GT.l.) WGRAD(l) = 0. 

WGRAD(NSP) = 0. 

GO TO 70 

65 LABEL! IM) = CHOKED 
70 CONI INUE 
FIRST = l 

I F ( AORB.GT. 1. ) FIRST = 2 
LAST = NSP Ml 

IF ( INTVL.GT.O) WRITE (6,1000) IM, I T 1 , ( WGRAD ( I ) , I =F I RST ,L AST ) 
80 IFIADRB.LE. 1.) RETURN 
WMB ( IM , 1 ) = WGRADt I ) 

WMB ( I M , 2 ) = WGRAD ( NS P ) 

WWCRMt IM, I) = WMB( IM, 1 I/WCR 
W WCRM( I M , 2 ) = WMB( IM, 2J/WCR 
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RETURN 

1000 FORMAT ( 5HK IM 
1010 FORMAT I73HK 
lVERTICAL LINE 
1020 FORMAT (92HK A 
150 ITERATIONS 
END 


. I 3, 10X , 5H ITl = , I 3/ I 2X, 10G 1 3 . A ) ) 

A VELOCITY GRADIENT SOLUTION CANNOT BE OBTAINED FOR 
IM =, 13 ) 

VELOCITY GRADIENT SOLUTION COULD NOT BE OBTAINED IN 
FOR VERTICAL LINE IM =,I3) 


BLOCK DATA 

COMMON /WWCRM/WWCRMI l 00, 2 ) , L ABEL ( 1 00 ) 

DATA L ABEL / l 00*6H 

END 


SUBROUTINE BLCD 
C 

C BLCC CALCULATES BLADE THETA COORDINATE AS A FUNCTION OF M 
C 

COMMON SRW, ITER, I END, LFR (2) ,NER(2> 

COMMON /INP/GAM, AR , T I P ,RHO I P , W TFL , OMEGA ,ORF , BE TA I , R E TAO , REDF AC , 

1 DENTOL » MB I * MBO, MM » NBB l » NBL, NRSP » MR ( 50 ) , RMSP ( 50 ) ,BESP(50) , 

2 BLDAT, AANDK, ERSOR , STRFN, SLCRD, INTVL.SURVL 
COMMON /C ALL ON/ ACTWT , ACT CMG , ACTL AM , MB IM I ,MEJ I PI ,MBOMI , MB OP I , MMM1 , 

1 HM1,HT,DTLR, OMLR, P I TCFi, C P , EXPON , TWW, CPT I P , TGROG , TB I , TBO, LAMBDA, 

2 TWL , ITMIN, ITMAX.NI P, IMS (2 ) ,BV{2) ,MV( 100) , IV 1 101 ) , I TV! 100,2) , 

3 TV! IOO,?),DTDMV(IOO,2),BETAV(100,2),MH( 100,2) »L)TUMH 1100,2), 

A B ET AH { 100,2) ,RMH( 1 CO , 2 ) , 6 EH { 1 00 , 2 ) , RM ( 100 ) , BE ( 1 00 ) , DB DM ( 1 00 ) , 

5 SAL ( 100 ), AAA ( 100) 

COMMON / G E 0 M I N / CHORD (2) , S TGR ( 2 ) , MLE ( 2 > , THLE 1 2 ) , RM I C2),RM0(2) , 

1 «I (?) ,R0(2) , BET II 2 ) .BETOI 2) .NSPII 2) .MSP l 50, 2) ,THSPI50,2) 

COMMON /BLCDCM/ EM { 50 , 2 ) , I N I T ( 2 ) 

COMMON /D2T0M 2/ 02TDM 2 ( 1 00 , 2 > 

ENTRY BLUM, THETA, DTCM, INF) 

INTEGER BLDAT , A ANDK, ER SOR , STRFN , SLCRD , SURVL , AA TEMP , SURF , F I RST, 

1 UPPER, SI, ST, SRW 

REAL K ,K AK ,LAMBDA,LMAX ,MH, MLE. MR , MSL .MSP, MV, MVIM I 
REAL M,MMLE,MSPMM,VMMSP 
SURF= I 
S I GN= 1. 

GO TO 10 

ENTRY BL21M, THETA, DTDM, INF) 

SI'RF= 2 
S I GN=— 1 . 

10 INF = 0 
IM = I 

CO 15 I =MB I , MBO 

15 IF ( ABS I MV ( I )-M) .LE.DMLR ) IM=I 
NSP- NS P I (SURF ) 

IF IINITISURFJ.EQ.13) GO TO 30 
ININSURF>= 13 

INITIAL CALCULATION OF FIRST AND LAST SPLINE POINTS ON BLADE 
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non 


c 

c 

c 


c 

c 

c 


A A = RETI ( SURF 1/57.295779 
A A - SINIAA) 

MSP ( L» SURF ) = R I ( SURF )* ( 1 .-SIGN*AA) 

B8 = SQRT ( 1. — AA**2) 

THSR( l * SURF ) = S IGN*RB*R I ( SURF) /RMI ( SURF) 

BET I (SURF ) = AA/BB/RMI ( SURF ) 

AA = BFTO(SURF)/57. 295779 
AA = SINIAA) 

MSP ( NSP » SURF ) = CHORD ( SURF ) -RO ( SURF ) * ( 1 . +S I GN* AA ) 

EB = SORT! l.-AA**2> 

THSPINSP, SURF ) = S TGR ( SURF ) + S I GN*BB»RO I SURF ) /RMO ( SURF ) 

BETH ( SURF ) = AA/BB/RMOI SURF ) 

CO 20 IA-1.NSP 

PSP ( IA, SURF ) = MSP( I A , SURF) +MLE C SURF ) 

20 THSP ( I A . SURF ) = THSP( I A , SURF ) +THLE ( SURF ) 

CALL S PLN22 ( MSP ( 1. SURF), THSP (1, SURF) ,BFTI (SURF ) , BE TO ( SURF ) , NS P , 

1 AAA, EM(1, SURF)) 

IF(DLDAT.LF.O) GO TO 30 
IF ( SURF. EG . I ) WR I TE ( 6, ICOO ) 

WRITE (6,1020) (MSP( I A. SURF), THSP (I A, SURF). AAA (I A). FM(I A, SURF), 

1 [ A= l , N S P ) 

BLACE COORDINATE C ALCUL A T ION 


30 KK = 2 

IF (M.GT.MSPl l, SURF) ) 
AT LEADING EDGE RADIUS 


GO TO 50 


AO 


M ML C = M-MLE(SURF) 

IF ( MMLE . LT.-DMLR ) GO TO 90 
MML E = AMAX110..MMLE) 

THETA- SQRT (MMLE* ( 2. *R I ( SURF J-MMLE ) ) *S IGN 

IF ( THETA . FQ.O . ) GO TO AO 

RHM= R I (SURF )-MML E 

CTDM= RMM/THETA/RM I ( SURF ) 

THETA = THE T A/RM I ( SURF ) 

C2TDM2( IM.SURF) = (-THET A-RMM*DTOM » / ( RM I ( SURF ) *THt TA ) **2 
THETA = THE T A+THLE (SURF ) 

RETURN 
INF- 1 

CTDM = 1 . E10 *S IGN 
THETA- THLE(SURF) 

C2TDM21 IM.SURF) = 0. 

RETURN 


ALONG SPLINE CURVE 

50 IF (M.LE.MSP(KK»SURF) ) GC TO 60 
IF (KK.GE.NSP) GO TO 70 
KK = KK+1 
GO TO 50 

60 S- MSP(KK»SURF)— MSPIKK— l.SURF) 
EMKM1- EM ( KK- L.SURF) 

E MK- E M ( KK , SURF ) 

MSPMM- MSPIKK, SURFJ-M 
M M MS P= M-MSP(KK-1, SURF ) 

THK- THSP(KK,SURF)/S 
THKM1- THSPIKK-l. SURF ) / S 
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c 

c 

c 


THE T A- EMKMl*MSPMM**3/6. /S + EMK*MMMSP»-»3/6./S «• 
l MMMSP + (THKMl-EMKMl*S/6. )*MSPMM 
CTQM= -EMKMl*MSPMM**2/2./S + EMK*MMMSP**2/2. /S * 
1 EMKM1 ) *S/6. 

C2TCM2 { I M, SURF ) = EMKM 1*MSPMM/S+EMK*MMMSP/S 
RETURN 


( THK-EMK * S /6 . )* 
THK-THKM1- ( EMK- 


AT TRAILING EDGE RADIUS 


70 


80 


CMM= CHORD ( SURFH-MLE ( SURF)-M 
IF (CMM.LT.-OMLR ) GO TO 90 
CMM = AMAXI (0. , CMM ) 

THETA= SORT I CMM* ( 2. *RG( SURF) -CMM) J*SIGN 
IF ( THETA. EQ.O. ) GO TC 80 
RMM = RO ( SURF )-CMM 


DTDM = -RMM/THETA/RMO ( SURF ) 

THETA = THETA/RMOI SURF ) 

C2TDM21 IM t SURF) = ( - THET A+RMM*D TDM ) / ( RMO l SURF I * THE T A ) **2 
THETA = THETA+STGR (SURF ) +THLF [ SURF ) 

RETURN 
I NF= 1 


DTDM = -1.EI0*SIGN 

THE TA = THEE ( SURF l+STGR(SURF J 

C2TDM2 ( I M , SURF ) » 0. 

RETURN 


ERRCR RETURN 


90 WRI TE(6, 1030) L ER ( 2 ) , V , SURF 
STOP 

1000 FORMAT ( 1 H l , 1 3 X » 33HBL AO E DATA AT INPUT SPLINE POINTS) 

1010 FORMAT! 1HL, 17X, 16HBLA0E SURFACE, 14) 

1020 FORMAT (7X , IHM, 10X, 5HTHET A , 10X , 10H0ER I VAT I VE , 5X , 1 0H2N0 DERIV. / 
l ( 4 G 1 5 • 5 > ) 

1030 FORMAT (IAHLBLCD CALL NO.,I3/33H M COORDINATE IS NOT HITHIN BLADE/ 
IAH M = ,G14.6, 10X.6HSURF =,G14.6) L Dt/ 

END 


FUNCTION IPF( IM, IT) 

COMMON /CALCON/ACTWT, ACTCMG, ACTLAM, MBIM1.MBIP1 ,MB0M1,MB0P1 

1 HM 1 , HT,nTLR,OMLR, PITCH, CP,EXPON,TWW,CPTIP, TGROG *TBI,T 0 O, LAMBDA , 

2 TWL » I TM IN , I TMAX ,N I P , I MS ( 2 ) , BV ( 2 ) , MV ( 100 ) , I V I 1 01 ) , I TVI 100 , 2 1 , 

3 IV( 100,2) ,DTDMV ( 100,2 ), BETAV! 100,2 ) , MH{ 100,2) ,DTDMH( 100,2) , 

4 BET AH ( 100,2), RMH( 1 00 , 2 ) , B EH ( 100, 2 ) , RM { 100 ) , BE I 1 00 ) , DBOM ( 1 00 ) . 

5 SAL( 100), AAAI 100) 

IPF = IV< I M ) + 1 T- I TV ( IM, 1 ) 

RETURN 

END 
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Subroutine CONTIN 


CONTIN calculates a new estimate for the initial value of W for equation (4). This 
is based on satisfying the continuity equation (7). If the input value of w (WTFL) is too 
large, there may be no solution of equation (4) satisfying equation (7). In this case, the 
choking weight flow will be found. 

An initial estimate of the velocity at the lower boundary is furnished by VELGRA, 
say W r The corresponding weight flow w^ is also calculated by VELGRA. CONTIN 
furnishes the next estimate W 2 , by linear interpolation or extrapolation from the origin 
(see fig. 16). Subsequent estimates are obtained by linear interpolation or extrapolation 
from the two previous estimates (see W 3 in fig. 16). This is essentially the method of 
false position (regula falsi). 



Figure 16. - Method used by subroutine CONTIN to determine relative hub velocity. 


If there is choked flow, so that a solution does not exist, information from three 
iterations is stored. This information is used to predict the next estimate of W on the 
lower boundary such that the weight flow will be a maximum. 

The input arguments for CONTIN are as follows: 

WA last value of W on lower boundary used in solving eq. (4) 

WTFL weight flow calculated by eq. (9) based on the input value of WA 

XND controls sequence of calculation in CONTIN; VELGRA sets END = 1 to 

indicate start of velocity- gradient solution for a new vertical mesh line 

I value of IM for vertical mesh line 

WT input weight flow, w 
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DELMAX maximum permitted change in estimated velocity WA per iteration 
TOLERC tolerance on velocity for calculating choking weight flow 

The output arguments for CONTEST are as follows: 


value of W on lower boundary to be used in solving eq. (4) 

used to control next iteration in CONTIN, and to indicate when a choked 
flow solution has been found 

The internal variables for CONTIN are as follows: 


DELTA predicted correction to WA 


NCALL number of iterations of CONTIN for a vertical mesh line 


SPEED 

WEIGHT 


array of values of WA from up to three previous iterations 
array of values of WTFL from up to three previous iterations 


SUBROUTINE CONTIN ( W A , WT FL , I NO , I , WT , DELMAX , TOL ERC ) 

C I MENS I UN SPEEOI 3) .WEIGHT! 3) 

NCALL = NCALL + L 

IF ( IND. NE.l. AND. NCALL. GT. 50) GO TO 400 
135 GO TO l 140, 150,210,270,370) , IND 
140 SPEED! 1 ) = WA 

WEIGHT! 1) = WTFL 
CELTA = WT/WTFL*WA-WA 

I F ( ABS ( DELTA ) . GT .DEL M AX ) DELTA = S I GN ( DELMAX , DEL TA ) 

IF ! wTFL.LT .0. ) DELTA ^ DELMAX 

WA = DELTA+WA 

IND = 2 

NCALL = 1 

RETURN 

150 IF ( (WTFL-WfcIGHT ( 1) I/IWA-SPEFD! 1) ) ) 180,180,160 
160 SPFFDJ2! = wA 

CELTA = ! WT-WTFL ) / ( WTFL-WE IGHT ( 1 ) ) * ( WA-SPEED 1 1 ) ) 
IF(ABS!DFLTA).GT. DELMAX) DELTA = S IGN I DE LMA X .DELTA ) 
WA = DELTA+WA 
166 SPEED! 1 ) = SPEFD ( 2 ) 

WEIGHT! 1) = WTFL 
RETURN 

170 WRITE 16,1000) I, WTFL, I 
IND = 6 
RETURN 
180 IND = 3 

IF (WTFL.GE.WT) GO TO 140 
IF (SPEED! I )-WA) 190 » 2C0 ,200 
190 SPEED! 2 ) = SPEED! 1 ) 

SPEED! 1) = 2.0*SPEED ( 1 )-WA 

S P E t D ( 3 ) = WA 

WEI GHT ( 2 ) = WE IGHT ( 1 ) 

WEIGHT! 3) = WTFL 
WA = SPEED! 1 ) 

RETURN 

200 SPEED! 2 ) » wA 



la: 

t mi = 

I* 
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SPEtO(3) = SPEED! 1) 

SPEED! 1 ) = 2 . 0*W A- SP E ED l 1) 

WEI GHT ( 2 ) = WTFL 
HEIGHT ( 3 ) = WEIGHTI1 ) 
toA = SPEED! 1 ) 

RETURN 

210 HEIGHT! 1) = WTFL 

IF (WTFL.GF.WT) GO TO 140 
IF (WEIGHT! 11-WEIGHT! 2)1 200,380,220 
220 WEI GHT { 3 ) = WEIGHT(2) 

WEI GHT ( 2 ) = W E I GHT ( 1 ) 

SPEED! 3 ) = SPEED! 2 ) 

SPEED! 2 ) = SPEED! 1 ) 

SPEED! 1 ) = 2.0*SPEED(2)-SPEED( 3) 

to A = SPEED! I ) 

RETURN 

230 IF (SPEED! 31-SPEED! ll-TCLERC) 170, 170, 240 

240 IND = 4 

245 IF (WE (GHT! 31-WEIGHT ( 1 ) ) 260,260,250 
250 toA = (SPEED! 1 ) + SPEED ( 2 ) 1/2.0 
RETURN 

260 toA = ( SPEED! 3) +SPEED! 2) 1 /2.0 
RETURN 

270 IFISPEFD! 31-SPEED! D-TCLERC) 170, 170, 280 

280 IF ( WT FL- WE I GHT ( 2 1 1 320,350,290 

290 IF (WA-SPEED! 2 I ) 310,300,300 

300 SPEED! 1 1 = SPEED ( 2 ) 

SPEED! 2 ) = WA 
toEIGHT! 1) = WEIGHT ( 2 ) 
toE I GHT ( 2 ) = WTFL 
GG TO 245 

310 SPEFDI3) = SPEED ( 2 ) 

SPEED! 2 ) = WA 
toE I GHT (3) = WE IGHT(2 1 
WEIGHT! 2) = WTFL 
GO TO 245 

320 IF ( WA-SPEED ( 2 ) ) 340,330,330 
330 WEIGHT! 3) = WTFL 
SPEED! 3 1 * to A 

GO TO 245 

340 WEIGHT! 1) = WTFL 
SPEED! I 1 = WA 

GO TO 245 
350 IND = 5 

IF ( WA-SPEED (2)1 380,360,360 

360 SPEED! 1 ) = SPEED ( 2 1 

WEIGHT! 1 ) = WEI GHT ( 2 1 

SPEED! 2 ) = ( SPEED ( 1 ) +SPEED ( 3 1 ) /2.0 

WA = SPEED ( 2 ) 

RETURN 
370 IND = 4 

WEI GHT ( 2 ) = WTFL 

WA ^ ( SPEED! 1 1 +SPEED! 2) ) /2.0 
RETURN 
380 IND = 5 

390 toEIGHT (3) = W£IGHT(2) 

SPEEDO) = SPEED ( 2 ) 

SPEED (2) = ( SPEED! 1) + SPEED! 3 > )/2. 

WA = SPEED! 2 ) 

RETURN 
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C NO SOLUTION FOUND IN 50 ITERATIONS 
400 IND = 7 
RETURN 

1000 FORMAT 1 43HLACTWT EXCEEDS CHOKING WEIGHT FLOW FUR IM =,13/ 
l 22HKCHOKING WEIGHT FLOW =,G15.6,9H FOR IM =,I3) 

END 


Subroutine INTGRL 

INTGRL calculator is the integral of a function passing through a given set of points. 
This subroutine is based on the spline curve. INTGRL solves a tridiagonal matrix equa- 
tion given in reference 9 to obtain the coefficients for the piecewise cubic polynomial 
function giving the spline fit curve. INTGRL is based on the end condition that the 
second derivative at either end point is one-half that at the next spline point. 

The input variables are as follows: 


X array of ordinates 

Y array of function values 

N integer number of X and Y values given 

The output variable is 

SUM array of values of integral of function, SUM(J) = Z ^ Y dX 

y X(l) 

If SRW = 17 in COMMON, input and output data for INTGRL are printed. This is 
useful in debugging. 


SUBROUTINE INTGRL (X,Y,N,SUM) 

C 

C INTGRL CALCULATES THE INTEGRAL OF A SPLINE CURVE PASSING THROUGH 
C A GIVEN SET OF POINTS 

C ENC CONDITION - SECOND DERIVATIVE AT EITHER END POINT IS ONE-HALF 
C THAT AT THE ADJACENT POINT 
C 

COMMON SRW 

COMMON / BOX/ G ( 50 ) , 50(50), EM(IOO) 

DIMENSION X ( N ) , Y(N), SUM(N) 

INTEGER SRW 
SB ( 1 ) = -.5 
G ( 1 ) = 0 
N0=N-l 

IFINO.LT.?) GO TO 20 
CO 10 1=2, Nu 
A = (XI I >-X( I-l))/6. 

C = (X( I + 1 )-X ( I ) )/6. 

W = 2 . * ( A +C ) -A* SB ( 1-1 ) 

SB (I) = C/W 

f = mi + n-Yim/ixd + n-xun-mn-Yii-m/ixm-xii-in 
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10 G ( I ) = ( F— A*G ( 1-1 ) )/W 

20 EMIN) = G ( N- l ) / ( 2 . +SB ( N- 1 ) ) 

CO 30 I =2 » N 
K = N+l-l 

30 EM l K ) = G(K)-SB(K)*EM(K+1) 

SUM ( 1 ) = 0.0 

CO jO K = 2,N 

50 SUM(K) = SUM ( K- l ) + l X ( K I — X C K- 1 ) ) * ( Y( K ) +Y IK- 1 ) ) /2 . 0- ( X I K) -X ( K 1)) 
l«*3*!EM(K)+EM(K-l) J/24.0 

I F ( SRW . EO . 17 ) WRITE! 6, 1000) N , ( X { I ) , Y ( I ) , SUM ( I ) , EM ( I ) , I = l , N ) 

RETURN 

1000 FORMAT ( l 7HK NO. OF PCINTS =I3/10X5HX 15X5HY 15X5HSUM 
l I3X10H2N0 OERIV./ ! 4E20.8) ) 

F NO 

The remaining subroutines are essentially the same as described in reference 1. 
The subroutines in TSONIC are not interchangeable with those in TANDEM or TURBLE , 
since there are differences in COMMON blocks and some changes in coding. However, 
the description of these subroutines in reference 1 still applies, with the exception of 
SPLINE and ROOT. In SPLINE the end condition has been changed so that the second 
derivative is the same at an end point as at the adjacent point. ROOT has been changed 
to find the root by the bisection method instead of by Newton's method. This is less 
efficient, but more foolproof. 


C 

c 

c 

c 

c 


SUBROUTINE SPLINE !X , Y . N .SLOPE , EM ) 


SPLINE CALCULATES FIRST AND SECOND DERIVATIVES AT SPL 1 ,N| ^ 

END CONDITION - SECOND DERIVATIVES ARE THE SAME AT END POINT AND 

ADJACENT POINT 


10 

20 


30 


COMMON SRW 

COMMON /BOX/ Gl 100) , SB ( ICO) 

C I MENS ION X ( N ) ,Y(N)»EM(N),SLOPE(N) 

INTEGER SRW 
SB! 1 ) = -1.0 
G ( 1 ) = 0 
N0=N— 1 

I F ( \0. L T . 2 ) GO TO 20 
CO 10 1=2. NO 
A = { X ( I ) — X ( 1-1) )/6. 

C = (X! 1 + 1 )-X( I ) )/6. 

W = 2. * { A+C ) -A*Sfl( I-l > 

S B { I ) = C / W 

F = I Y ( 1 + 1 )-Y( I ) ) / !X ( l + l )— X I I) )- I Y ( I )-YI I- 1 ) ) / (X l I ) -X (1-1 ) ) 
G ( I ) = ! F — A » G ( I-l) )/W 

EM!N) = GIN— l)/< l.+SBIN-l) ) 


CO 30 1=2, N 
K = N+l-I 

SLO^Et I ) G = K (XtU-X(2| l |/6!*(2.»EMt 1)*FM(21 l/.XUJ-XUM 


AO SLORE(I) = ! XI I ) -X ( I -l ) ) /6 . » l 2. »EM ( I ) +EM( I-l ) ) M Y( l ) Y(I 1))/ 

1 IF! SRW ! EQ • 1 3 ) 1 WRITE (6,1000) N, (XII), Y(I) .SLOPE (I),EM(l),I = l.N) 

RETURN 


r> o n o 


1000 FORMAT (2X.15HNU. OF POINTS 
1 2HEM/ ( 4G20 . 8 ) ) 

END 


= , I 3/1 OX, 1HX, 19X, IHY, 19X, 5H SLOPE , 1SX , 


SUBROUTINE SPLN22 (X , Y , Y IP , YNP ,N, SLOPE , EM J 

CAt - CuL fl TES FIRST AND SECOND DERIVATIVES AT SPL IN C POINTS 
END CONDITION - DERIVATIVES SPECIFIED AT END POINTS 

COMMON SRW 

COMMON /BOX/ G 1 100 ) , SB ( 1 00 ) 

CIMENSION X(N),Y(N),EN(N),SLOPE(N) 

INTEGER SRW 
SB ( 1 ) = .5 

F = l Y ( 2)-Y( 1 ) )/ (X(2 )-X( 1 ) l-YIP 
G( 1 ) = F*3./ ( X(2 ) — X ( 1 ) ) 

N0=N-1 

IF ( NO.LT. 2 ) GU TO 20 
CO 10 1=2, NO 
A = (X ( I )-X( I — 1 1 )/6. 

C = (X( I+l)-X( I) )/6. 

W = 2. * ( A+C ) -A»SB ( 1-1 ) 

SB ( I ) = C/W 

20 F = YNP-IY(N)-YIN-l) )/(X(N)-X(N-l) ) 

W = (X(N)-X(N-i) )/6.*(2.-SB(N-l) ) 

EMIN) = (F-(X(N)-X(N-l ) ) *G(N-1 )/6. )/W 
CO 30 1=2, N 
K = N+l-I 

30 EM(K) = G(K)-SB(K)*EM(K+1) 

co°4o ( il2lN tX(U_X(2 ,,/6,M2,fEM,i,tEM[2 i> + [ Y ^)-Ym)/(xi2i-xmi 
^0^ SCOPE ( I )^= ( ( X ( I )— X ( I — 1 ) ) /6 . * ( 2 .*EM ( I ) -*-EM ( I - 1 J J + J Y ( I ) — Y ( I - 1 ) ) / 

RETURN* 60 ’ 18> WRIT£ (6 ’ 1COO) ^*^XI I),Y(I), SLOPE (I),EM(I),I = l f N) 

1000 12SEMJuG20!ln N °* ° F P0If<TS = » I3/l0X ' 1HX «l 9 ^lHY,19X,5HSL0PE,15X, 
END 


SUBROUTINE SPLINT I X , Y ,N ,Z , M AX , YI NT , DYDX ) 

C 

C SPLINT CALCULATES INTERPOLATED POINTS AND DERIVATIVES 
C FOR A SPLINE CURVE 

C ENC CONDITION - SECOND DERIVATIVE AT CITHER END POINT IS ONE-HALF 
C THAT AT THE ADJACENT POINT 
C 

COMMON SRW 

COMMON /BOX/ Gl 100), SB! 100) 

CIMENSION X(N) ,Y(N),Z (MAX) ,Y INTI MAX) .DYDX (MAX) 

CIMENSION EM { 100 ) 


EQUIVALENCE (SB, EM) 

INTEGER SRW 

I F ( MAX • LE . 0 ) RETURN 

III = SRW 

SB ( 1 ) = - • 5 

G( I) = 0 

N0=N-1 

I F ( NO* LT • 2 ) GO TO 20 
CO 10 I =2 « NO 
A = (X ( I >-X( 1-1 ) )/6. 

C = (XU + 1)-X( I ) )/6. 

W = 2* * ( A + C ) -A*SB { 1-1 ) 

SB ( I ) = C/W 

f = (Y( i+i )-Y( m/tx c i+n-xc 1 1 )-cy( n-Yc i-i i ) / tx( i )-xc i-in 
10 G ( I ) = (F-A*G( 1-1) )/W 
20 EMIN) = GIN-l)/(2.+SB(N-L> ) 

CO 30 I = 2 * N 
K = N+l-I 

30 EM ( K ) = G(K)-SB(K) *EM(K+1) 

DO l AO 1 = 1 t MAX 
K = 2 

IF(Z(I)-Xll)) 70,60,90 
60 Y INT ( I ) =Y ( 1 ) 

SK = X { K ) — X ( K— 1 ) 

GO TO 130 

70 IF(ZU).GE.(l.l*X(l)-.l*X<2> )) GO TO 120 
WRITE (6,1000) ZII) 

SRW - 16 
GO TO 120 
80 K = N 

IF(Z(I> .LE. { 1. 1*X(N>-. 1*X(N-1) )) GO TO 120 
WRITE (6,1000) Z(I) 

SRW = 16 
GO TO 120 

90 I F ( Z ( I )-X(K) ) 120,100,110 

100 YINTU ) = Y ( K ) 

SK = X(K>-X{K-1) 

GO TO 130 
110 K=K+ 1 

IF(K-N) 90,90,80 
120 CONTINUE 

SK = X(K)—X(K— 1) 

Y I NT ( I ) = EM(K-1 )* (X ( K)-Z( I ) )**3/6. /SK +EM ( K ) * ( Z (I)-X ( K- 1 ) ) ** 3/6. 

1 / SK+ ( Y ( K ) /SK -EM(K)*SK / 6. ) * ( Z ( I ) -X ( K- l ) ) + ( Y ( K- 1 ) / SK -EM(K-1 ) 

2 *SK/6.)*U(K)-Z( I) ) 

130 CYDX ( I )=-EM(K-l )*(X(K )-Z ( I))**2/2.0/SK +EM ( K ) * ( X ( K- 1 ) -Z (I) ) **2/2. 

1 / SK+ ( Y { K ) -Y ( K- 1 ) ) /SK - ( EM ( K ) -EM ( K- 1 ) )*SK/6. 

I AO CONTINUE 

MXA = MAXO ( N , M AX ) 

I F ( SRW. EQ. 16) WRITE! 6, 1010) N, MAX , ( X ( I) , Y (I) , Z (I) , Y I NT (I ) , DYDX ( I ) , 
1 1=1, MXA) 

SRW = III 
RETURN 

1000 FORMAT ( 5AH SPLINT USED FOR EXTRAPOLATION. EXTRAPOLATED VALUE = * 
1G1A.6) 

1010 FORMAT ( 2 X , 2 1HN0. OF POINTS GIVEN =,I3,30H f NO. OF INTERPOLATED PO 
1 I NTS =, I3/10X, 1HX , 19X , 1HY, 16X, 1 IHX- 1 NTERPOL . . 9X, 11HY- INTERPOL. , 
28X, 1AHDY0X- INTERPOL. / ( 5E 20. 8) > 

END 
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SUBROUTINE MHURI Z ( MV , I TV , BL , MB I ,MBO, I TO, HT , 0 IX R , KOOF f J, MH, DTDMH, 
1MRTS) 

MHORIZ CALCULATES M COORDINATES OF INTERSECTIONS OF ALL HORIZONTAL 
ME SF LINES WITH A BLADF SURFACE 
KODE = 0 FOR UPPER BLADE SURFACE 
KODF = l FOR LOWER BLADE SURFACE 

COMMON SRW, ITER, I END, LERI2 J , NFR( 2) 

DIMENSION MV ( 100 ) , IT V 1 100) ,MH( 100 ), DTDMH l 100) 

INTEGER BLDAT, AANDK, ERSOR, STKFN, SLCRD , SUR VL , AA TEMP , SURF , F I R ST , 

1 UPPER, SI, ST, SRW 

REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MV IMl 
REAL MV I M 
EXTERNAL BL 

IF (MBI.GE.MBO) RETURN 
I M= MB I 
10 ITIND* 0 

20 IF ( I T V ( IM+lJ-ITVl IMJ-ITINDJ 30,40,50 
30 J= J+l 

T I = FLOAT ( I TV ( IM+1 )- ITC-IT IND + KODE ) *HT 
IT I N0= IT I N 0 - 1 
MV I M = MV ( IM ) 

IF (MRTS.EQ.I) MVIM = M V IM+ ( MV ( l M+ 1 ) -MVI M ) / l 000. 

CALL ROOT (MVIM,MV(IM+1) ,TI,BL,DTLP. , MH ( J ) , DTDMH ( J ) ) 

CO TO 20 
40 I M= IM+I 
MRTS = 0 

IF (IM.EQ.MBO) RETURN 
CO TO 10 
50 J= J+l 

T I = FLOAT! ITV( IMI-ITO+IT INC+KOOE ) *HT 
I T I s iD= ITINU+1 
MVI K = MVI IM) 

IF (MRTS.EQ.I) MVIM = MV IM + ( MV I I M+ 1 ) -MV I M ) / 1 000. 

CALL ROOT (MVIM ,MV( I M + 1 ) , T I , BL , DTLR , MH< J ) , D TUMH { J ) ) 

GO TO 20 
END 


SUBROUTINE DEN ST Y ( RHOw , RHO , V EL , T WLMR , CPT I P , E XPUN , RHO I P , G AM , AR , T I P ) 
C 

C DENSTY CALCULATES DENSITY AND VELOCITY FROM THE WEIGHT FLOW PARAMETER 
C DENSITY TIMES VFLOCITY 
C 

COMMON SRW, I TER, IEND,LFR (2) ,NER( 2) 

VEL = RHOW/RHO 
IF (VEL.NE.O. ) GO TO 10 
RHO = RHO IP 
RETURN 

10 T T I 0 = l.-( VEL**2+TWLMR ) /CPTIP 
IF ( TTIP.LT.O. ) GO TO 30 
TEM^ = TT I P* * ( EX PON- 1 . ) 

RHOr = RHOIP*TEMP*TT IP 

R HOW.P_= ; - VEL** 2/GAM* RHO I P/ AR *T EMP/ T I P+RHOT 
IF (RHOWP. LE.O. ) GO TO 30 
VELNEW = VEL+ ( RHOW-RHCT *VEL ) /RHOWP 
IF(A8S(VELNE W- VEL) /V ELNEW.LT. .0001) GO TO 20 
VEL = VELNEW 
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GO TO 10 
20 VEL = VELNEw 
RHO = RHOW/VEL 
RETURN 

30 T GROG = 2.*GAM*AR/(GAM+1,) 

VEL = SQRT(TGROG*TIP*( l.-TWLMR/CPTIP) ) 

RHO = RHOIPM l.-(VFL**2+TWLMR)/CPTIP)**EXP0N 
RWMURW = RHGW/RHO/VEL 
NER(l) = NER(i)+l 

WR I T F ( 6 , 1000 ) LER(1> , NER ( 1 ) , RWMQft W 

I F ( NFR ( 1) .FQ.50) STOP 

RETURN 

1000 FORMAT ( 16HL0ENSTY CALL NC.,I3/9H NER { 1 ) =,I3/10H RHC*W IS ,F7.A, 
1 3 AH TIMES THE MAXIMUM VALUE FOR RHO*W) 

ENO 


SUBROUTINE ROOT ( A , ft, Y f FU NCT , TOL ERY , X , DFX ) 

ROOT FINDS A ROOT FOR (FUNCT MINUS Y) IN THE INTERVAL (A,8) 

COMMON SRWtITERt I END » LER ( 2 ) »NER(2) 

INTEGER SRW 

IF (SRW.EQ.21) WR I TE { 6 f I COO ) A , B , Y , TOLERY 
XI = A 

CALL FUNCTU1, FXltOFX, INF) 

IFISRW.EQ.21 ) WRITE(6. 1010) X 1 , FX 1 , DFX , INF 
X2 = B 

10 CO 30 1=1,20 
X = (X1+X2 )/2. 

CALL FLJNCTIX.FX, DFX, INF) 

IFC SRW.EQ.21 ) WRITEI6, 1010) X f FX,DFX,INF 
I F ( (FX1-Y)*<FX-Y) .GT.C. ) GO TO 20 
X2 = X 
GO TO 30 
20 XI = X 
FX1 = FX 
30 CONTINUE 

I F ( ABS ( Y-FX ) .LT. TOLERY ) RETURN 
WRITE! 6, 1020) LER(2),A,B,Y 
STOP 

1000 FORMAT ( 3 2 H 1 INPUT ARGUMENTS FOR ROOT — A =G 13 . 5 , 3X , 3HB =,G13.S, 

1 3X.3HY =,G1 3. 5, 3X, 8HTGLERY = ,G1 3. 5/16X, IHX, 17X , 2HFX , 1 5X, 3HDFX , 

2 1 OX , 3H I NF ) 

1010 F0RMAT(8X,G16.5»2G18.5, 16) 

1020 FORMAT ( 1AHLR00T CALL NC.,I3/37H ROOT HAS FAILED TO OBTAIN VALID RO 
1CT/AH A = ,G1A.6, 10X, 3HB = , G 1 A . 6 , 1 0 X , 3HY =,G1A.6) 

END 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, June 27, 1969, 
720-03-00-66-22. 


APPENDIX A 


DERIVATION OF VELOCITY-GRADIENT EQUATION 

The velocity-gradient equation is an expression of the force equation. By a balance 
of force in the ^-direction, the following equation can be obtained: 

0W_ 1 d ( rV 0) _ 1 d(rW 0 +a>r 2 ) 

00 W dt W dt * A1 

The time derivative indicates the change in the quantity for a moving particle as a func- 
tion of time. Equation (Al) is a special case of equation (BIO) of reference 11. We 
make use of the following relations (see fig. 1): 


W 0 = W sin p 


W m = W cos (3 

W r = W m sin a 

W - cos a 
z m 


— = w 

dt r 


r 


d0 

dt 


= W 0 


dm = w 

dt 


m 



dt dS 

+ w 0W 
dt y r00 3m 
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We can perform the indicated differentiation of the right side of equation (Al) by 
using the preceding relations , and then solve for 3W/ 9 0 (which appears on both sides 
of the equation) to obtain 

= Si + sin a tan A + r tan 0 55 + 2o>r «l2-2L (A2) 

30 \cos /3 dS / 9m cos /3 

We desire now to evaluate dp/dS in terms of first and second derivatives of 0 
with respect to m along streamlines. The following relations hold for streamlines: 

tan 0 = — (A3) 

dm 

r9u 

tan /3 = — = - SHH. (A4) 

W m iH 

30 

Equation (A4) is obtained by using equations (2) and (3). Also along streamlines we have 

= cos 0 (A 5) 

dS 


— = sin a 
dm 


(A6) 


Now differentiate equation (A3) and use equations (A5) and (A6) to obtain 


d/3 = r cos 3 p d2g i sin a sin A cos2 A 


dS 


dm 


(A7) 


2 2 

Along the surface of the blade d 0/dm can be easily calculated since 0 is given 

O P 

explicitly as a function of m. However, in the passage d 0/dm is given indirectly 

p P 

by the stream function. Hence, we will need an expression for d 0/dm in terms of the 
partial derivatives of the stream function. First, from equations (A3) and (A4) we have 

3u 

d0 _ 3m (A8) 

dm 3u 

90 
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By differentiating equation (A8) we obtain 


d 2 6> = _3_ /_d0_\ + _3_ /_d0\ dd_ 
dm 2 3m \dm/ 56 \dm / dm 


2 3u _3u_ 8 2 u _ / 3u\ 2 3 2 u 

30 3m 80 3m \30/ 3 2 

dm 




(A9) 


Finally, by using the fact (from eqs. (A3) and (A8)) that 

r cos 0 _ sin 0 
8u 8u 

30 3m 


we can obtain 


r 


2 





3u 

30 


(A10) 


By using equations (A7) and (A10), equation (A2) can be put in the following form: 


— = AW+B (All) 

30 


where 


2 2 d 2 0 9 

A = r cos /3 + sin a tan 0(1 + cos^ 0) 


is used on blade surface, 


drn 


(A12a) 


A = sin 2 0 


a2 

8 u 

30 3m 


8u 

30 


^2 
3 u 


A 

30 2 


iH. / 3u\ 2 
3m \3m) 


3m 2 — 

30 


+ sin a tan 0(1 + cos z 0) 


(A12b) 


90 


is used at interior points, and 


B = r tan /3 — + 2u>r sm - (A13) 

3m cos /3 

Equations (All) to (A13) are in the form used in the program. 
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APPENDIX B 


DEFINING REDUCED WEIGHT FLOW PROBLEM 

Since the final solution obtained by the velocity-gradient equation depends on the 
stream-function solution obtained with a reduced weight flow, it is important to establish 
the conditions which will give the most suitable stream-function solution. To accomplish 
this, the streamlines at the reduced weight flow should correspond in shape as closely 
as possible to those for the actual weight flow. Some of the factors affecting this corre- 
spondence are discussed in the following paragraphs. 

The first condition is that equation (1) should not be changed. This condition is 
satisfied if the ratio w/w is not changed. Therefore, is reduced in the same ratio 
as w, for the reduced weight flow solution. With this condition satisfied, there would 
be no change at all in the stream function, if the flow were incompressible. With com- 
pressible flow, the stream -function solution will change since the coefficients in equa- 
tion (1) are functions of the density p, which in turn is a function of the relative velocity 
W. The relative velocity naturally will be lower if the weight flow is reduced. 

Another consideration is the boundary conditions at the upstream and downstream 
boundaries, AH and DE (see fig. 4). We could use the same boundary condition as for 
the full weight flow. However, this is not the best approximation. The reason for 
this is that the input information is the mean value of 0 le and p at BG and CF, 
respectively, instead of AH and DE. And we want to obtain a streamline pattern sat- 
isfying the input condition, but at a reduced weight flow. So the entire calculation of 
^in and ^out for the stream -function solution is based on the reduced weight flow 
condition. Notice, also, that the calculation of /3^ n and j3 Qut requires a value of 
prerotation A, which in turn depends on co and w. Hence, a value of A based on 
the reduced weight flow is also calculated. The method of calculating A and p. and 
^out is described in appendix B of reference 1. 

We can summarize the quantities which must be calculated based on the reduced 
weight flow. They are w, «, A, p and 0 out . These quantities, based on the reduced 
weight flow, are used for calculating the reduced weight flow solution. Values of w, w, 
and A based on the actual weight flow are also calculated for later use in the velocity- 
gradient equations. 
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